Utilities: Locate and square/circle intersection area

.github/workflows/ci.yml

@@ -120,6 +120,7 @@ jobs:
       - name: Set Environment Variables
         run: |
           echo "PELE_PHYSICS_HOME=${{github.workspace}}/PelePhysics-${{matrix.comp}}" >> $GITHUB_ENV
+          echo "UTILITIES_WORKING_DIRECTORY=${{github.workspace}}/PelePhysics-${{matrix.comp}}/Testing/Exec/UtilityUnitTest" >> $GITHUB_ENV
           echo "TRANSPORT_WORKING_DIRECTORY=${{github.workspace}}/PelePhysics-${{matrix.comp}}/Testing/Exec/TranEval" >> $GITHUB_ENV
           echo "EOS_WORKING_DIRECTORY=${{github.workspace}}/PelePhysics-${{matrix.comp}}/Testing/Exec/EosEval" >> $GITHUB_ENV
           echo "REACT_WORKING_DIRECTORY=${{github.workspace}}/PelePhysics-${{matrix.comp}}/Testing/Exec/ReactEval" >> $GITHUB_ENV
@@ -174,6 +175,24 @@ jobs:
             cd ${{github.workspace}}/PelePhysics-${{matrix.comp}}/Support/ceptr
             poetry install
           fi
+      - name: Unit Test Utilities
+        working-directory: ${{env.UTILITIES_WORKING_DIRECTORY}}
+        run: |
+          echo "::add-matcher::${{github.workspace}}/PelePhysics-${{matrix.comp}}/.github/problem-matchers/gcc.json"
+          if [ "${{matrix.comp}}" == 'hip' ]; then source /etc/profile.d/rocm.sh; fi;
+          if [ "${{matrix.comp}}" == 'sycl' ]; then source /opt/intel/oneapi/setvars.sh || true; fi;
+          ccache -z
+          make -j TINY_PROFILE=TRUE USE_CCACHE=TRUE ${{matrix.amrex_build_args}}; \
+          if [ "${{matrix.comp}}" == 'gnu' ] || [ "${{matrix.comp}}" == 'llvm' ]; then \
+            ./PelePhysicsUnitTest2d.${{matrix.comp}}.TPROF.ex;
+          fi;
+          make realclean;
+          if [ $? -ne 0 ]; then exit 1; fi; \
+      - name: Utilities ccache report
+        working-directory: ${{env.UTILITIES_WORKING_DIRECTORY}}
+        run: |
+          ccache -s
+          du -hs ${HOME}/.cache/ccache
       - name: Test Transport
         working-directory: ${{env.TRANSPORT_WORKING_DIRECTORY}}
         run: |

Docs/sphinx/Utility.rst

@@ -294,7 +294,9 @@ The user must ensure that the correct number of grow cells is present in the Fab
 
 Utilities
 =============================
-The utilities namespace contains unit conversions, which are particularly useful for PeleLM(eX) users working with mixed unit systems. The following aliases, defined in a file header, enable straightforward conversions between MKS and CGS units for use in equation of state (``eos``) function calls:
+The ``utilities`` namespace contains some functions that may be useful to users when writing problem-specific source code.
+
+This includes unit conversions, which are particularly useful for PeleLM(eX) users working with mixed unit systems. The following aliases, defined in a file header, enable straightforward conversions between MKS and CGS units for use in equation of state (``eos``) function calls:
 
 ::
 
@@ -314,3 +316,6 @@ For example, users can call ``eos.PYT2R`` as follows:
    // Calculate density in CGS units, then convert to MKS
    eos.PYT2R(m2c::P(P_mean), massfrac, Temp, rho_cgs);
    amrex::Real rho = c2m::Rho(rho_cgs); // Convert eos density to MKS
+
+The ``utilities`` namespace also contains some other useful functions: ``locate`` to find the nearby indices of a value in an array for interpolation
+and ``rectangle_circle_intersection_area`` to find the area of the intersections between rectangles and circles.

Source/Utility/BlackBoxFunction/Table.H

@@ -5,6 +5,7 @@
 #include <iostream>
 #include <fstream>
 #include "BlackBoxFunction.H"
+#include "Utilities.H"
 
 namespace pele::physics {
 
@@ -469,48 +470,6 @@ private:
   const TabulatedFunctionData* tf_data;
   static constexpr int maxdim = TableDim != 0 ? TableDim : MAXD_TABLE;
 
-  // -----------------------------------------------------------
-  // Search for the closest index in an array to a given value
-  // using the bisection technique.
-  // INPUTS/OUTPUTS:
-  // xtable(0:n-1) => array to search in (ascending order)
-  // n             => array size
-  // x             => x location
-  // idxlo (return)=> output st. xtable(idxlo) <= x < xtable(idxlo+1)
-  // -----------------------------------------------------------
-  AMREX_GPU_HOST_DEVICE
-  AMREX_FORCE_INLINE
-  int locate(const amrex::Real* xtable, const int n, const amrex::Real x)
-  {
-    int idxlo = 0;
-    // If x is out of bounds, return boundary index
-    if (x >= xtable[n - 2]) {
-      idxlo = n - 2;
-      return idxlo;
-    }
-    if (x <= xtable[0]) {
-      idxlo = 0;
-      return idxlo;
-    }
-
-    // Do the bisection
-    int idxhi = n - 2;
-    bool notdone = true;
-    while (notdone) {
-      if (idxhi - idxlo <= 1) {
-        notdone = false;
-      } else {
-        const int idxmid = (idxhi + idxlo) / 2;
-        if (x >= xtable[idxmid]) {
-          idxlo = idxmid;
-        } else {
-          idxhi = idxmid;
-        }
-      }
-    }
-    return idxlo;
-  }
-
   // Get quantities for use in interpolation
   AMREX_GPU_HOST_DEVICE
   AMREX_FORCE_INLINE
@@ -530,18 +489,14 @@ private:
       auto* beg = &tf_data->grids[start];
       int idx = static_cast<int>(
         std::lower_bound(beg + 1, beg + dimLen - 1, interpdata[dim]) - beg - 1);
-      idx = std::min(
-        std::max(idx, 0),
-        dimLen - 2); // don't let idx be negative or out of bounds
 #else
-      // home rolled
-      // int idx = 1;
-      // while (tf_data->grids[start+idx] < interpdata[dim] && idx < dimLen)
-      //    {idx+=1;};
-      // idx-=1;
-      int idx = locate(&tf_data->grids[start], dimLen, interpdata[dim]);
+      int idx = 0;
+      utilities::locate(&tf_data->grids[start], dimLen, interpdata[dim], idx);
 #endif
 
+      idx = std::min(
+        std::max(idx, 0),
+        dimLen - 2); // don't let idx be negative or out of bounds
       indices[dim] = idx;
       idx += start;
       amrex::Real dx_inv =

Source/Utility/Utilities/Utilities.H

@@ -4,6 +4,124 @@
 #include "UnitConversions.H"
 namespace pele::physics::utilities {
 
+// -----------------------------------------------------------
+// Compute the overlap area between an arbitrary rectangle and circle
+// code modified from:
+// https://stackoverflow.com/questions/622287/area-of-intersection-between-circle-and-rectangle
+// INPUTS:
+// x0, x1        => low and high x values of rectangle
+// y0, y1        => low and high y values of rectangle
+// cx, cy        => x and y values of circle center
+// r             => radius of circle
+// OUTPUT:
+// area of intersection in unit/coordinate systemm of inputs
+// -----------------------------------------------------------
+
+// helper function
+AMREX_GPU_HOST_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+section_area(
+  const amrex::Real x0,
+  const amrex::Real x1,
+  const amrex::Real h,
+  const amrex::Real r)
+{
+  const amrex::Real sec = (h < r) ? std::sqrt(r * r - h * h) : 0;
+  const amrex::Real arg1 = amrex::max(-sec, amrex::min(sec, x1));
+  const amrex::Real arg2 = amrex::max(-sec, amrex::min(sec, x0));
+  const amrex::Real term1 =
+    0.5 * (std::sqrt(1.0 - arg1 * arg1 / (r * r)) * arg1 * r +
+           r * r * std::asin(arg1 / r) - 2.0 * h * arg1);
+  const amrex::Real term2 =
+    0.5 * (std::sqrt(1.0 - arg2 * arg2 / (r * r)) * arg2 * r +
+           r * r * std::asin(arg2 / r) - 2.0 * h * arg2);
+  return term1 - term2;
+}
+
+// Area of a square and circle intersection
+AMREX_GPU_HOST_DEVICE
+AMREX_FORCE_INLINE
+amrex::Real
+rectangle_circle_intersection_area(
+  const amrex::Real x0,
+  const amrex::Real x1,
+  const amrex::Real y0,
+  const amrex::Real y1,
+  const amrex::Real cx,
+  const amrex::Real cy,
+  const amrex::Real r)
+{
+  // We assume coordinates of the box are passed in sorted order x0 = xlo, x1 =
+  // xhi, etc
+  AMREX_ASSERT(x1 > x0);
+  AMREX_ASSERT(y1 > y0);
+  AMREX_ASSERT(r >= 0.0);
+
+  // Convert to coordinates centered on circle
+  const amrex::Real x0new = x0 - cx;
+  const amrex::Real x1new = x1 - cx;
+  const amrex::Real y0new = y0 - cy;
+  const amrex::Real y1new = y1 - cy;
+
+  // Return 0 if fully outside the square that circumscribes the circle
+  if (x0new >= r || x1new <= -r || y0new >= r || y1new <= -r) {
+    return 0.0;
+  }
+
+  // Compute the area
+  const amrex::Real sign = copysign(1.0, y0new) * copysign(1.0, y1new);
+  const amrex::Real y0mod = amrex::min(std::abs(y0new), std::abs(y1new));
+  const amrex::Real y1mod = amrex::max(std::abs(y0new), std::abs(y1new));
+  const amrex::Real area0 = section_area(x0new, x1new, 0.0, r);
+  const amrex::Real area1 = area0 - section_area(x0new, x1new, y0mod, r);
+  const amrex::Real area2 = area0 - section_area(x0new, x1new, y1mod, r);
+
+  return area2 - sign * area1;
+}
+
+// -----------------------------------------------------------
+// Search for the closest index in an array to a given value
+// using the bisection technique.
+// INPUTS/OUTPUTS:
+// xtable(0:n-1) => array to search in (ascending order)
+// n             => array size
+// x             => x location
+// idxlo        <=> output st. xtable(idxlo) <= x < xtable(idxlo+1)
+// -----------------------------------------------------------
+AMREX_GPU_HOST_DEVICE
+AMREX_FORCE_INLINE
+void
+locate(const amrex::Real* xtable, const int n, const amrex::Real& x, int& idxlo)
+{
+  // If x is out of bounds, return boundary index
+  if (x >= xtable[n - 1]) {
+    idxlo = n - 1;
+    return;
+  }
+  if (x <= xtable[0]) {
+    idxlo = 0;
+    return;
+  }
+
+  // Do the bisection
+  idxlo = 0;
+  int idxhi = n - 1;
+  bool notdone = true;
+  while (notdone) {
+    if (idxhi - idxlo <= 1) {
+      notdone = false;
+    } else {
+      const int idxmid = (idxhi + idxlo) / 2;
+      if (x >= xtable[idxmid]) {
+        idxlo = idxmid;
+      } else {
+        idxhi = idxmid;
+      }
+    }
+  }
+}
+
 } // namespace pele::physics::utilities
 
-#endif
+#endif

Testing/Exec/UtilityUnitTest/GNUmakefile

@@ -0,0 +1,78 @@
+# AMReX
+DIM          = 2
+PRECISION    = DOUBLE
+PROFILE      = FALSE
+TINY_PROFILE = FALSE
+VERBOSE      = FALSE
+DEBUG        = FALSE
+
+# Compiler
+COMP       = gnu
+FCOMP      = gfortran
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+USE_CUDA   = FALSE
+USE_HIP    = FALSE
+
+PELE_PHYSICS_HOME ?= ../../..
+AMREX_HOME        ?= $(abspath $(lastword $(PELE_PHYSICS_HOME)/Submodules/amrex))
+UTILITY_HOME      ?= $(PELE_PHYSICS_HOME)/Source/Utility
+
+EBASE = PelePhysicsUnitTest
+BL_NO_FORT = TRUE
+
+##
+
+Bpack   := ./Make.package
+Blocs   := .
+
+PP_SRC_HOME = $(PELE_PHYSICS_HOME)/Source
+EXTERN_CORE       += $(PP_SRC_HOME)
+INCLUDE_LOCATIONS += $(PP_SRC_HOME)
+VPATH_LOCATIONS   += $(PP_SRC_HOME)
+Bpack             += $(PP_SRC_HOME)/Make.package
+Blocs             += $(PP_SRC_HOME)
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+INCLUDE_LOCATIONS += $(UTILITY_HOME)
+VPATH_LOCATIONS   += $(UTILITY_HOME)
+#include $(UTILITY_HOME)/Diagnostics/Make.package
+#include $(UTILITY_HOME)/Filter/Make.package
+#include $(UTILITY_HOME)/PMF/Make.package # Needs mechanism.H but shouldn't
+#include $(UTILITY_HOME)/PltFileManager/Make.package
+#include $(UTILITY_HOME)/TurbInflow/Make.package
+#include $(UTILITY_HOME)/BlackBoxFunction/Make.package
+include $(UTILITY_HOME)/Utilities/Make.package
+
+Pdirs   := Base AmrCore Boundary
+Bpack	+= $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
+Blocs	+= $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir))
+include $(Bpack)
+INCLUDE_LOCATIONS += $(Blocs)
+VPATH_LOCATIONS   += $(Blocs)
+
+PP_SRC_HOME = $(PELE_PHYSICS_HOME)/Source
+EXTERN_CORE       += $(PP_SRC_HOME)
+INCLUDE_LOCATIONS += $(PP_SRC_HOME)
+VPATH_LOCATIONS   += $(PP_SRC_HOME)
+Bpack             += $(PP_SRC_HOME)/Make.package
+Blocs             += $(PP_SRC_HOME)
+
+# job_info support
+CEXE_sources += AMReX_buildInfo.cpp
+CEXE_headers += $(AMREX_HOME)/Tools/C_scripts/buildInfo.H
+INCLUDE_LOCATIONS +=  $(AMREX_HOME)/Tools/C_scripts
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
+
+clean::
+	$(SILENT) $(RM) AMReX_buildInfo.cpp
+
+AMReX_buildInfo.cpp:
+	$(AMREX_HOME)/Tools/C_scripts/makebuildinfo_C.py \
+          --amrex_home "$(AMREX_HOME)" \
+          --COMP "$(COMP)" --COMP_VERSION "$(COMP_VERSION)" \
+          --FCOMP "$(FCOMP)" --FCOMP_VERSION "$(FCOMP_VERSION)" \
+          --MODULES "$(MNAMES)" \
+          --GIT "$(PELE_PHYSICS_HOME) $(AMREX_HOME) $(SUNDIALS_HOME)"

Testing/Exec/UtilityUnitTest/Make.package

@@ -0,0 +1 @@
+CEXE_sources += main.cpp