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