clean up linear interpolation

.github/workflows/model_parser.yml

@@ -0,0 +1,34 @@
+name: model parser
+
+on: [pull_request]
+jobs:
+  model_parser:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+
+      - name: Get submodules
+        run: |
+          git submodule update --init
+          cd external/Microphysics
+          git fetch; git checkout development
+          cd ../amrex
+          git fetch; git checkout development
+          cd ../..
+
+      - name: Install dependencies
+        run: |
+          sudo apt-get update -y -qq
+          sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0
+
+      - name: Compile model_parser test
+        run: |
+          cd Util/model_parser/test
+          make -j 4
+
+      - name: Run model_parser test
+        run: |
+          cd Util/model_parser/test
+          ./Castro3d.gnu.ex

Util/model_parser/model_parser.H

@@ -64,6 +64,10 @@ namespace model_string
 }
 
 
+///
+/// return the index into the model coordinate, loc, such that
+/// model::profile(model_index).r(loc) < r < model::profile(model_index).r(loc+1)
+///
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
 int
 locate(const Real r, const int model_index) {
@@ -74,7 +78,7 @@ locate(const Real r, const int model_index) {
        loc = 0;
 
     } else if (r > model::profile(model_index).r(model::npts-2)) {
-       loc = model::npts-1;
+       loc = model::npts-2;
 
     } else {
 
@@ -94,18 +98,14 @@ locate(const Real r, const int model_index) {
         loc = ihi;
     }
 
-    return loc;
+    return loc-1;
 }
 
 
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
 Real
 interpolate(const Real r, const int var_index, const int model_index=0) {
 
-    // find the value of model_state component var_index at point r
-    // using linear interpolation.  Eventually, we can do something
-    // fancier here.
-
     int id = locate(r, model_index);
 
     Real slope;
@@ -119,12 +119,11 @@ interpolate(const Real r, const int var_index, const int model_index=0) {
        interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);
 
        // safety check to make sure interp lies within the bounding points
-       Real minvar = amrex::min(model::profile(model_index).state(id+1, var_index),
+       Real minvar = std::min(model::profile(model_index).state(id+1, var_index),
                                 model::profile(model_index).state(id, var_index));
-       Real maxvar = amrex::max(model::profile(model_index).state(id+1, var_index),
+       Real maxvar = std::max(model::profile(model_index).state(id+1, var_index),
                                 model::profile(model_index).state(id, var_index));
-       interp = amrex::max(interp, minvar);
-       interp = amrex::min(interp, maxvar);
+       interp = std::clamp(interp, minvar, maxvar);
 
     } else if (id == model::npts-1) {
 
@@ -135,12 +134,11 @@ interpolate(const Real r, const int var_index, const int model_index=0) {
 
 
        // safety check to make sure interp lies within the bounding points
-       Real minvar = amrex::min(model::profile(model_index).state(id-1, var_index),
+       Real minvar = std::min(model::profile(model_index).state(id-1, var_index),
                                 model::profile(model_index).state(id, var_index));
-       Real maxvar = amrex::max(model::profile(model_index).state(id-1, var_index),
+       Real maxvar = std::max(model::profile(model_index).state(id-1, var_index),
                                 model::profile(model_index).state(id, var_index));
-       interp = amrex::max(interp, minvar);
-       interp = amrex::min(interp, maxvar);
+       interp = std::clamp(interp, minvar, maxvar);
 
     } else {
 
@@ -166,9 +164,11 @@ interpolate(const Real r, const int var_index, const int model_index=0) {
 
 }
 
-// Subsample the interpolation to get an averaged profile. For this we need to know the
-// 3D coordinate (relative to the model center) and cell size.
 
+///
+/// Subsample the interpolation to get an averaged profile. For this we need to know the
+/// 3D coordinate (relative to the model center) and cell size.
+///
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real interpolate_3d (const Real* loc, const Real* dx, int var_index, int nsub = 1, int model_index = 0)
 {

Util/model_parser/test/GNUmakefile

@@ -0,0 +1,33 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 3
+
+COMP	   = gnu
+
+USE_MPI    = FALSE
+USE_OMP    = FALSE
+
+USE_ALL_CASTRO = FALSE
+USE_AMR_CORE = FALSE
+
+# define the location of the CASTRO top directory
+CASTRO_HOME  ?= ../../..
+
+USE_MODEL_PARSER = TRUE
+NUM_MODELS := 2
+
+# This sets the EOS directory in Castro/EOS
+EOS_DIR     := helmholtz
+
+# This sets the network directory in Castro/Networks
+NETWORK_DIR := subch_simple
+
+EXTERN_SEARCH += .
+
+Bpack   := ./Make.package
+Blocs   := . .. $(CASTRO_HOME)/Source/problems
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Util/model_parser/test/Make.package

@@ -0,0 +1,4 @@
+CEXE_sources += main.cpp
+CEXE_sources += extern_parameters.cpp
+CEXE_headers += extern_parameters.H
+

Util/model_parser/test/README.md

@@ -0,0 +1,4 @@
+# Unit test for model_parser
+
+This simply reads in an initial model and does some checks to make
+sure the locate and interpolation routines work as expected.

Util/model_parser/test/main.cpp

@@ -0,0 +1,41 @@
+#include <iostream>
+
+#include <model_parser.H>
+
+
+int main(int argc, char *argv[]) {
+
+    amrex::Initialize(argc, argv);
+
+    // initialize the external runtime parameters in C++
+
+    init_extern_parameters();
+
+    // now initialize the C++ Microphysics
+#ifdef REACTIONS
+    network_init();
+#endif
+
+    Real small_temp = 1.e-200;
+    Real small_dens = 1.e-200;
+    eos_init(small_temp, small_dens);
+
+    std::string model = "sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.N.10.00km";
+    read_model_file(model);
+
+    Real r{3.89e7};
+
+    std::cout << "testing locate" << std::endl;
+
+    auto idx = locate(r, 0);
+    AMREX_ALWAYS_ASSERT(r >= model::profile(0).r(idx) &&
+                        r <= model::profile(0).r(idx+1));
+
+    std::cout << "testing interpolate" << std::endl;
+
+    // density is monotonically decreasing
+    auto dens = interpolate(r, model::idens);
+    AMREX_ALWAYS_ASSERT(dens <= model::profile(0).state(idx, model::idens) &&
+                        dens >= model::profile(0).state(idx+1, model::idens));
+
+}