fix race condition in slopez for 4th order case near physical boundaries

update the treatment to match slopex and slopey cases which recompute slopes intead of relying on neighboring slopes which may not have been computed yet on a GPU
AMReX-Astro · Dec 6, 2023 · fea9749 · fea9749
1 parent 80c3e1a
commit fea9749
Showing 1 changed file with 13 additions and 3 deletions.
diff --git a/Source/MaestroSlopes.cpp b/Source/MaestroSlopes.cpp
@@ -537,12 +537,17 @@ void Maestro::Slopez(const Box& bx, Array4<Real> const s,
 
                     } else if (k == klo + 1) {
                         // Recalculate the slope at lo(2)+1 using the revised dzl
+                        Real del = -16.0 / 15.0 * s(i, j, k - 2, n) +
+                                   0.5 * s(i, j, k - 1, n) +
+                                   2.0 / 3.0 * s(i, j, k, n) -
+                                   0.1 * s(i, j, k + 1, n);
                         dmin = 2.0 * (s(i, j, k - 1, n) - s(i, j, k - 2, n));
                         dpls = 2.0 * (s(i, j, k, n) - s(i, j, k - 1, n));
                         Real slim = amrex::min(amrex::Math::abs(dpls),
                                                amrex::Math::abs(dmin));
                         slim = dpls * dmin > 0.0 ? slim : 0.0;
-                        dzl = slz(i, j, k - 1, n);
+                        Real sflag = amrex::Math::copysign(1.0, del);
+                        dzl = sflag * amrex::min(slim, amrex::Math::abs(del));
                         ds = 4.0 / 3.0 * dcen - (dzr + dzl) / 6.0;
                         slz(i, j, k, n) =
                             dflag * amrex::min(amrex::Math::abs(ds), dlim);
@@ -568,12 +573,17 @@ void Maestro::Slopez(const Box& bx, Array4<Real> const s,
 
                     } else if (k == khi - 1) {
                         // Recalculate the slope at lo(3)+1 using the revised dzr
+                        Real del = -(-16.0 / 15.0 * s(i, j, k + 2, n) +
+                                     0.5 * s(i, j, k + 1, n) +
+                                     2.0 / 3.0 * s(i, j, k, n) -
+                                     0.1 * s(i, j, k - 1, n));
                         dmin = 2.0 * (s(i, j, k + 1, n) - s(i, j, k, n));
                         dpls = 2.0 * (s(i, j, k + 2, n) - s(i, j, k + 1, n));
                         Real slim = amrex::min(amrex::Math::abs(dpls),
                                                amrex::Math::abs(dmin));
                         slim = dpls * dmin > 0.0 ? slim : 0.0;
-                        dzr = slz(i, j, k + 1, n);
+                        Real sflag = amrex::Math::copysign(1.0, del);
+                        dzr = sflag * amrex::min(slim, amrex::Math::abs(del));
                         ds = 4.0 / 3.0 * dcen - (dzl + dzr) / 6.0;
                         slz(i, j, k, n) =
                             dflag * amrex::min(amrex::Math::abs(ds), dlim);
@@ -582,4 +592,4 @@ void Maestro::Slopez(const Box& bx, Array4<Real> const s,
             });
     }
 }
-#endif
+#endif