diff --git a/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp b/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp index d2d7272f91d..8e344ec572b 100644 --- a/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp +++ b/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp @@ -14,6 +14,8 @@ #include #include +#include +#include #include #include #include @@ -357,6 +359,54 @@ namespace amrex::Abort("SetNeigh: dim must be 0, 1 or 2"); } + + /** + * \brief Get the address of the value of arr in the neighbor (i_n, j_n) on + * the plane with normal 'dim' (same indexing convention as GetNeigh), for + * atomic updates of the neighbor's value. + * + * \param[in] arr data to be accessed + * \param[in] i, j, k the indices of the "center" cell + * \param[in] i_n the offset of the neighbor in the first direction + * \param[in] j_n the offset of the neighbor in the second direction + * \param[in] dim normal direction to the plane in consideration (0 for x, 1 for y, 2 for z) + */ + template + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + constexpr + T* + GetNeighPtr(const amrex::Array4& arr, + const int i, const int j, const int k, + const int i_n, const int j_n, const int dim){ + + if(dim == 0){ + return &arr(i, j + i_n, k + j_n); + } + #ifdef WARPX_DIM_XZ + else if(dim == 1 || (dim == 2)){ + return &arr(i + i_n, j + j_n, k); + } + #elif defined(WARPX_DIM_3D) + else if(dim == 1){ + return &arr(i + i_n, j, k + j_n); + } + else if(dim == 2){ + return &arr(i + i_n, j + j_n, k); + } + #else + else if(dim == 1){ + amrex::Abort("GetNeighPtr: Only implemented in 2D3V and 3D3V"); + } + else if(dim == 2){ + return &arr(i + i_n, j + j_n, k); + } + #endif + + amrex::Abort("GetNeighPtr: dim must be 0, 1 or 2"); + + return nullptr; + } + #ifdef AMREX_USE_EB #ifndef WARPX_DIM_RZ @@ -676,17 +726,33 @@ WarpX::ComputeOneWayExtensions () // has given away already some area, so we use Sz_red rather than Sz. // If no face is available we don't do anything and we will need to use the // multi-face extensions. - const int flag_neigh = - ::GetNeigh(flag_info_face, i, j, k, i_n, j_n, idim); - if (::GetNeigh(S_mod, i, j, k, i_n, j_n, idim) > S_ext - && (flag_neigh == FaceInfo::available - || flag_neigh == FaceInfo::intruded) - && flag_ext_face(i, j, k)) { - - ::SetNeigh(S_mod, - ::GetNeigh(S_mod, i, j, k, i_n, j_n, idim) - S_ext, - i, j, k, i_n, j_n, idim); - + // The area is taken with an atomic test-and-subtract: on GPU + // several faces can try to borrow from the same intruded face + // concurrently, and a plain read-test-write lets the intruded + // face give the same area away more than once. + // + // *Atomically* decrement `S_mod` of the neighboring cell by + // `S_ext` under the condition that this is possible (i.e. that + // S_mod-S_ext is positive, and the cell is marked as available + // for borrowing). If this indeed updated `S_mod`, it returns + // `true` for `borrowed`. For the syntax, see + // https://amrex-codes.github.io/amrex/doxygen/namespaceamrex_1_1Gpu_1_1Atomic.html + const bool borrowed = amrex::Gpu::Atomic::If( + ::GetNeighPtr(S_mod, i, j, k, i_n, j_n, idim), // address to atomically update + S_ext, // value to combine + amrex::Minus(), // operation to perform when combining the value + // condition: callable that gets called with rem=S_mod-S_ext. + // The `flag_ext_face` test also stops this loop after the + // first successful borrow, as it is cleared just below. + [=] (amrex::Real rem) { + const int flag_neigh = ::GetNeigh(flag_info_face, i, j, k, i_n, j_n, idim); + return rem > amrex::Real(0.) + && (flag_neigh == FaceInfo::available + || flag_neigh == FaceInfo::intruded) + && flag_ext_face(i, j, k); + }); + + if (borrowed) { // Insert the index of the face info borrowing_inds[ps] = ps; // Store the information about the intruded face in the dataset of the @@ -705,6 +771,17 @@ WarpX::ComputeOneWayExtensions () } } } + // The counting pass reserved one slot for this face, but the atomic + // test-and-subtract above fails if a concurrently extended face drained + // the intruded face in the meantime. The face is then still flagged, and + // has to report that it borrowed nothing: the solver reads + // `borrowing_size` entries starting at `*borrowing_inds_pointer` (see + // EvolveBCartesianECT), which were never filled in. The face itself is + // left to the eight-ways extension. + if (flag_ext_face(i, j, k)) { + borrowing_size(i, j, k) = 0; + borrowing_inds_pointer(i, j, k) = nullptr; + } } }, amrex::Scan::Type::exclusive); } @@ -854,10 +931,36 @@ WarpX::ComputeEightWaysExtensions () if(denom >= S_ext){ S_mod(i, j, k) = S(i, j, k); int count = 0; + // The extension is all-or-nothing: a face that got only some of its + // patches would not reach its stable area, and the area it did take + // would be lost to the faces that lent it, since the ECT update of an + // intruded face assumes that the area it lent is accounted for by the + // face that borrowed it (see EvolveBCartesianECT). + bool all_borrowed = true; for (int i_n = -1; i_n < 2; i_n++) { for (int j_n = -1; j_n < 2; j_n++) { if(local_avail(i_n + 1, j_n + 1) != 0_rt){ + if (count == nborrow) { + // The borrowing pass found more available neighbors + // than the counting pass reserved slots for + all_borrowed = false; + continue; + } const amrex::Real patch = S_ext * ::GetNeigh(S, i, j, k, i_n, j_n, idim) / denom; + // Atomic test-and-subtract, for the same reason as in + // ComputeOneWayExtensions: an intruded face shared by + // concurrently extended faces must not give the same + // area away more than once + const bool borrowed = amrex::Gpu::Atomic::If( + ::GetNeighPtr(S_mod, i, j, k, i_n, j_n, idim), + patch, amrex::Minus(), + [=] (amrex::Real rem) { + return rem > amrex::Real(0.); + }); + if (!borrowed) { + all_borrowed = false; + continue; + } borrowing_inds[ps + count] = ps + count; FaceInfoBox::addConnectedNeighbor(i_n, j_n, ps + count, borrowing_neigh_faces); @@ -868,14 +971,39 @@ WarpX::ComputeEightWaysExtensions () i, j, k, i_n, j_n, idim); S_mod(i, j, k) += patch; - ::SetNeigh(S_mod, - ::GetNeigh(S_mod, i, j, k, i_n, j_n, idim) - patch, - i, j, k, i_n, j_n, idim); count +=1; } } } - flag_ext_face(i, j, k) = false; + if (!all_borrowed) { + // Give the area of the successful patches back, and leave the face + // flagged so that it is stabilized by the BCK correction instead + for (int n = 0; n < count; n++) { + auto const vec = + FaceInfoBox::uint8_to_inds(borrowing_neigh_faces[ps + n]); + amrex::Gpu::Atomic::AddNoRet( + ::GetNeighPtr(S_mod, i, j, k, vec(0), vec(1), idim), + borrowing_area[ps + n]); + } + count = 0; + S_mod(i, j, k) = S(i, j, k); + } + // The recorded size has to match the entries actually written, since + // the solver reads `borrowing_size` entries starting at + // `*borrowing_inds_pointer` (see EvolveBCartesianECT) + borrowing_size(i, j, k) = count; + if (count == 0) { + borrowing_inds_pointer(i, j, k) = nullptr; + } else { + flag_ext_face(i, j, k) = false; + } + } + else { + // The area available shrank between the counting and the borrowing + // pass: the face cannot be extended after all, and has to report + // that it borrowed nothing + borrowing_size(i, j, k) = 0; + borrowing_inds_pointer(i, j, k) = nullptr; } } }, amrex::Scan::Type::exclusive);