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 9072 : BCondBase<Field>::BCondBase(unsigned int face)
24 9072 : : face_m(face)
25 9072 : , 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 672 : 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 672 : unsigned int face = this->face_m;
41 672 : unsigned d = face / 2;
42 672 : if (field.getCommunicator().size() > 1) {
43 672 : const Layout_t& layout = field.getLayout();
44 672 : const auto& lDomains = layout.getHostLocalDomains();
45 672 : const auto& domain = layout.getDomain();
46 672 : int myRank = field.getCommunicator().rank();
47 :
48 672 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
49 768 : || (lDomains[myRank][d].min() == domain[d].min());
50 :
51 672 : if (!isBoundary) {
52 0 : return;
53 : }
54 672 : }
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 672 : typename Field::view_type& view = field.getView();
61 672 : 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 672 : if (nghost > 1) {
67 0 : throw IpplException("ExtrapolateFace::apply", "nghost > 1 not supported");
68 : }
69 :
70 672 : 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 672 : if (face & 1) {
76 336 : src = view.extent(d) - 2;
77 336 : dest = src + 1;
78 : } else {
79 336 : src = 1;
80 336 : 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 3584 : for (unsigned i = 0; i < Dim; i++) {
87 2912 : begin[i] = nghost;
88 2912 : end[i] = view.extent(i) - nghost;
89 : }
90 672 : begin[d] = src;
91 672 : end[d] = src + 1;
92 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
93 672 : ippl::parallel_for(
94 672 : "Assign extrapolate BC", createRangePolicy<Dim, exec_space>(begin, end),
95 4492672 : KOKKOS_CLASS_LAMBDA(index_array_type & args) {
96 : // to avoid ambiguity with the member function
97 : using ippl::apply;
98 :
99 2245664 : T value = apply(view, args);
100 :
101 2245664 : args[d] = dest;
102 :
103 2245664 : apply(view, args) = slope_m * value + offset_m;
104 : });
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 1680 : void PeriodicFace<Field>::findBCNeighbors(Field& field) {
139 1680 : auto& comm = field.getCommunicator();
140 : // For cell centering only face neighbors are needed
141 1680 : unsigned int face = this->face_m;
142 1680 : unsigned int d = face / 2;
143 1680 : const int nghost = field.getNghost();
144 1680 : int myRank = comm.rank();
145 1680 : const Layout_t& layout = field.getLayout();
146 1680 : const auto& lDomains = layout.getHostLocalDomains();
147 1680 : const auto& domain = layout.getDomain();
148 :
149 16240 : for (auto& neighbors : faceNeighbors_m) {
150 14560 : neighbors.clear();
151 : }
152 :
153 1680 : if (lDomains[myRank][d].length() < domain[d].length()) {
154 : // Only along this dimension we need communication.
155 :
156 480 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
157 720 : || (lDomains[myRank][d].min() == domain[d].min());
158 :
159 480 : if (isBoundary) {
160 : // this face is on mesh/physical boundary
161 : // get my local box
162 480 : auto& nd = lDomains[myRank];
163 :
164 : // grow the box by nghost cells in dimension d of face
165 480 : auto gnd = nd.grow(nghost, d);
166 :
167 : int offset;
168 480 : if (face & 1) {
169 : // upper face
170 240 : offset = -domain[d].length();
171 : } else {
172 : // lower face
173 240 : offset = domain[d].length();
174 : }
175 : // shift by offset
176 480 : gnd[d] = gnd[d] + offset;
177 :
178 : // Now, we are ready to intersect
179 1440 : for (int rank = 0; rank < comm.size(); ++rank) {
180 960 : if (rank == myRank) {
181 480 : continue;
182 : }
183 :
184 960 : if (gnd.touches(lDomains[rank])) {
185 240 : faceNeighbors_m[face].push_back(rank);
186 : }
187 : }
188 : }
189 : }
190 1680 : }
191 :
192 : template <typename Field>
193 672 : void PeriodicFace<Field>::apply(Field& field) {
194 672 : auto& comm = field.getCommunicator();
195 672 : unsigned int face = this->face_m;
196 672 : unsigned int d = face / 2;
197 672 : typename Field::view_type& view = field.getView();
198 672 : const Layout_t& layout = field.getLayout();
199 672 : const int nghost = field.getNghost();
200 672 : int myRank = comm.rank();
201 672 : const auto& lDomains = layout.getHostLocalDomains();
202 672 : const auto& domain = layout.getDomain();
203 :
204 : // We have to put tag here so that the matchtag inside
205 : // the if is proper.
206 672 : int tag = comm.next_tag(mpi::tag::BC_PARALLEL_PERIODIC, mpi::tag::BC_CYCLE);
207 :
208 672 : if (lDomains[myRank][d].length() < domain[d].length()) {
209 : // Only along this dimension we need communication.
210 :
211 192 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
212 288 : || (lDomains[myRank][d].min() == domain[d].min());
213 :
214 192 : if (isBoundary) {
215 : // this face is on mesh/physical boundary
216 : // get my local box
217 192 : auto& nd = lDomains[myRank];
218 :
219 : int offset, offsetRecv, matchtag;
220 192 : if (face & 1) {
221 : // upper face
222 96 : offset = -domain[d].length();
223 96 : offsetRecv = nghost;
224 96 : matchtag = comm.preceding_tag(mpi::tag::BC_PARALLEL_PERIODIC);
225 : } else {
226 : // lower face
227 96 : offset = domain[d].length();
228 96 : offsetRecv = -nghost;
229 96 : matchtag = comm.following_tag(mpi::tag::BC_PARALLEL_PERIODIC);
230 : }
231 :
232 192 : 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 192 : 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 192 : HaloCells_t& halo = field.getHalo();
241 192 : std::vector<range_t> rangeNeighbors;
242 :
243 288 : for (size_t i = 0; i < neighbors.size(); ++i) {
244 96 : int rank = neighbors[i];
245 :
246 96 : auto ndNeighbor = lDomains[rank];
247 96 : ndNeighbor[d] = ndNeighbor[d] - offset;
248 :
249 96 : NDIndex<Dim> gndNeighbor = ndNeighbor.grow(nghost, d);
250 :
251 96 : NDIndex<Dim> overlap = gndNeighbor.intersect(nd);
252 :
253 : range_t range;
254 :
255 432 : for (size_t j = 0; j < Dim; ++j) {
256 336 : range.lo[j] = overlap[j].first() - nd[j].first() + nghost;
257 336 : range.hi[j] = overlap[j].last() - nd[j].first() + nghost + 1;
258 : }
259 :
260 96 : rangeNeighbors.push_back(range);
261 :
262 : detail::size_type nSends;
263 96 : halo.pack(range, view, haloData_m, nSends);
264 :
265 96 : buffer_type buf = comm.template getBuffer<memory_space, T>(nSends);
266 :
267 96 : comm.isend(rank, tag, haloData_m, *buf, requests[i], nSends);
268 96 : buf->resetWritePos();
269 : }
270 :
271 288 : for (size_t i = 0; i < neighbors.size(); ++i) {
272 96 : int rank = neighbors[i];
273 :
274 96 : range_t range = rangeNeighbors[i];
275 :
276 96 : range.lo[d] = range.lo[d] + offsetRecv;
277 96 : range.hi[d] = range.hi[d] + offsetRecv;
278 :
279 96 : detail::size_type nRecvs = range.size();
280 :
281 96 : buffer_type buf = comm.template getBuffer<memory_space, T>(nRecvs);
282 96 : comm.recv(rank, matchtag, haloData_m, *buf, nRecvs * sizeof(T), nRecvs);
283 96 : buf->resetReadPos();
284 :
285 : using assign_t = typename HaloCells_t::assign;
286 96 : halo.template unpack<assign_t>(range, view, haloData_m);
287 : }
288 192 : if (!requests.empty()) {
289 96 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
290 : }
291 192 : comm.freeAllBuffers();
292 192 : }
293 : // For all other processors do nothing
294 : } else {
295 480 : if (d >= Dim) {
296 0 : throw IpplException("PeriodicFace::apply", "face number wrong");
297 : }
298 :
299 480 : 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 2720 : for (size_t i = 0; i < Dim; ++i) {
309 2240 : end[i] = view.extent(i) - nghost;
310 2240 : begin[i] = nghost;
311 : }
312 480 : begin[d] = 0;
313 480 : end[d] = nghost;
314 :
315 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
316 480 : ippl::parallel_for(
317 480 : "Assign periodic field BC", createRangePolicy<Dim, exec_space>(begin, end),
318 4253632 : 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 672 : }
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
|
