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 ${{env.NPROCS}} TINY_PROFILE=TRUE USE_CCACHE=TRUE ${{matrix.amrex_build_args}}; \
+          if [ "${{matrix.comp}}" == 'gnu' ] || [ "${{matrix.comp}}" == 'llvm' ]; then \
+            ./Pele2d.${{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,44 @@
+# AMReX
+DIM        = 2
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+VERBOSE    = FALSE
+DEBUG      = FALSE
+
+# Compiler
+COMP       = gnu
+FCOMP      = gfortran
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+USE_CUDA   = FALSE
+USE_HIP    = FALSE
+
+# PelePhysics
+FUEGO_GAS  = TRUE
+TINY_PROFILE = FALSE
+
+# define the location of the PELE_PHYSICS top directory
+PELE_PHYSICS_HOME    ?= ../../..
+
+# this flag activates the subcycling mode in the D/Cvode routines
+DEFINES  += -DMOD_REACTOR
+
+# Activates use of SUNDIALS
+ifeq ($(USE_CUDA), TRUE)
+  PELE_USE_KLU = FALSE
+else
+  ifeq ($(USE_HIP), TRUE)
+    PELE_USE_KLU = FALSE
+  else
+    PELE_USE_KLU = FALSE
+  endif
+endif
+
+Eos_Model       = GammaLaw
+Chemistry_Model = Null
+Transport_Model = Constant
+
+Bpack   := ./Make.package
+Blocs   := .
+
+include $(PELE_PHYSICS_HOME)/Testing/Exec/Make.PelePhysics

Testing/Exec/UtilityUnitTest/Make.package

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

Testing/Exec/UtilityUnitTest/main.cpp

@@ -0,0 +1,62 @@
+#include <AMReX.H>
+#include <AMReX_Print.H>
+#include "Utilities.H"
+
+namespace m2c = pele::physics::utilities::mks2cgs;
+namespace c2m = pele::physics::utilities::cgs2mks;
+
+int
+main(int argc, char* argv[])
+{
+  amrex::Initialize(argc, argv);
+
+  // Test Unit Conversions
+  AMREX_ALWAYS_ASSERT(m2c::Length(1.0) == 100.0);
+  AMREX_ALWAYS_ASSERT(c2m::Length(1.0, 2) == 1.0e-4);
+  AMREX_ALWAYS_ASSERT(m2c::Mass(c2m::Mass(1.0)) == 1.0);
+  AMREX_ALWAYS_ASSERT(m2c::Mu(1.0) == m2c::Mass(1.0) / m2c::Length(1.0));
+  amrex::Print() << "Unit Conversion Tests Passed" << std::endl;
+
+  // Test locate function
+  amrex::Array<amrex::Real, 5> sorted_data{0.0, 0.2, 0.5, 0.9, 1.0};
+  int length = sorted_data.size();
+  int idxlo = -1;
+  amrex::Real x = -50.0;
+  pele::physics::utilities::locate(sorted_data.data(), length, x, idxlo);
+  AMREX_ALWAYS_ASSERT(idxlo == 0);
+  x = 1.0;
+  pele::physics::utilities::locate(sorted_data.data(), length, x, idxlo);
+  AMREX_ALWAYS_ASSERT(idxlo == length - 1);
+  x = 0.25;
+  pele::physics::utilities::locate(sorted_data.data(), length, x, idxlo);
+  AMREX_ALWAYS_ASSERT(idxlo == 1);
+  amrex::Print() << "Locate Tests Passed" << std::endl;
+
+  // Test rectangle circle intersection area
+  // Quarter circle intersection
+  amrex::Real x0 = 0.0;
+  amrex::Real x1 = 1.0;
+  amrex::Real y0 = 0.0;
+  amrex::Real y1 = 1.0;
+  amrex::Real cx = 1.0;
+  amrex::Real cy = 1.0;
+  amrex::Real r = 1.0;
+  amrex::Real A1 = pele::physics::utilities::rectangle_circle_intersection_area(
+    x0, x1, y0, y1, cx, cy, r);
+  AMREX_ALWAYS_ASSERT(std::abs(A1 - 0.25 * M_PI) < 1e-14);
+  // No intersection
+  cx = 3.0;
+  cy = 3.0;
+  amrex::Real A2 = pele::physics::utilities::rectangle_circle_intersection_area(
+    x0, x1, y0, y1, cx, cy, r);
+  AMREX_ALWAYS_ASSERT(std::abs(A2) < 1e-14);
+  // Square covered by circle
+  r = 20.0;
+  amrex::Real A3 = pele::physics::utilities::rectangle_circle_intersection_area(
+    x0, x1, y0, y1, cx, cy, r);
+  AMREX_ALWAYS_ASSERT(std::abs(A3 - 1.0) < 1e-14);
+  amrex::Print() << "Intersection Area Tests Passed" << std::endl;
+
+  amrex::Finalize();
+  return 0;
+}