Merge branch 'stable' into dev

AFD-Illinois · Aug 28, 2024 · e5c777b · e5c777b
2 parents a7bd87b + da55d42
commit e5c777b
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 10 deletions.
diff --git a/kharma/reductions/reductions.hpp b/kharma/reductions/reductions.hpp
@@ -51,7 +51,9 @@ namespace Reductions {
 std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages);
 
 /**
- * Perform a reduction using operation 'op' over a given domain
+ * Perform a reduction using operation 'op' over a given domain.
+ * Note startx/stopx are in *embedding* coordinates e.g. Kerr-Schild, not MKS/FMKS
+ * 
  * This should be used for all 2D shell sums not around the EH:
  * Just set equal min/max, 2D slices are detected
  */

diff --git a/kharma/reductions/reductions_impl.hpp b/kharma/reductions/reductions_impl.hpp
@@ -119,9 +119,9 @@ T Reductions::CheckOnAll(MeshData<Real> *md, int channel)
 }
 
 #define INSIDE (x[1] > startx1 && x[2] > startx2 && x[3] > startx3) && \
-                (trivial1 ? x[1] < startx1 + G.Dxc<1>(i) : x[1] < stopx1) && \
-                (trivial2 ? x[2] < startx2 + G.Dxc<2>(j) : x[2] < stopx2) && \
-                (trivial3 ? x[3] < startx3 + G.Dxc<3>(k) : x[3] < stopx3)
+                (trivial1 ? xin[1] < startx1 : x[1] < stopx1) && \
+                (trivial2 ? xin[2] < startx2 : x[2] < stopx2) && \
+                (trivial3 ? xin[3] < startx3 : x[3] < stopx3)
 
 // TODO additionally template on return type to avoid counting flags with Reals
 template<Reductions::Var var, typename T>
@@ -173,11 +173,13 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
         pmb0->par_reduce("domain_sum", block.s, block.e, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
             KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
                 const auto& G = U.GetCoords(b);
-                GReal x[4];
+                GReal x[GR_DIM], xin[GR_DIM];
                 G.coord_embed(k, j, i, Loci::center, x);
+                if (trivial1 || trivial2 || trivial3)
+                    G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
                 if(INSIDE) {
                     local_result += reduction_var<var>(REDUCE_FUNCTION_CALL) *
-                        (!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
+                        ((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
                 }
             }
         , sum_reducer);
@@ -189,11 +191,13 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
         pmb0->par_reduce("domain_max", block.s, block.e, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
             KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
                 const auto& G = U.GetCoords(b);
-                GReal x[4];
+                GReal x[GR_DIM], xin[GR_DIM];
                 G.coord_embed(k, j, i, Loci::center, x);
+                if (trivial1 || trivial2 || trivial3)
+                    G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
                 if(INSIDE) {
                     const Real val = reduction_var<var>(REDUCE_FUNCTION_CALL) *
-                        (!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
+                        ((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
                     if (val > local_result) local_result = val;
                 }
             }
@@ -206,11 +210,13 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
         pmb0->par_reduce("domain_min", block.s, block.e, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
             KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
                 const auto& G = U.GetCoords(b);
-                GReal x[4];
+                GReal x[GR_DIM], xin[GR_DIM];
                 G.coord_embed(k, j, i, Loci::center, x);
+                if (trivial1 || trivial2 || trivial3)
+                    G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
                 if(INSIDE) {
                     const Real val = reduction_var<var>(REDUCE_FUNCTION_CALL) *
-                        (!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
+                        ((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
                     if (val < local_result) local_result = val;
                 }
             }