AMReX-Astro · zingale · Dec 26, 2023 · Sep 17, 2023 · Sep 26, 2023 · Oct 18, 2023
diff --git a/.github/workflows/model_parser.yml b/.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
diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out
@@ -1,5 +1,5 @@
 #   COLUMN 1                        2                        3                        4                        5                        6                        7                        8                        9                       10                       11                       12                       13                       14                       15                       16                       17                       18                       19                       20                       21                       22                       23
 #   TIMESTEP                     TIME                     MASS                     XMOM                     YMOM                     ZMOM              ANG. MOM. X              ANG. MOM. Y              ANG. MOM. Z              KIN. ENERGY              INT. ENERGY               GAS ENERGY             GRAV. ENERGY             TOTAL ENERGY     CENTER OF MASS X-LOC     CENTER OF MASS Y-LOC     CENTER OF MASS Z-LOC     CENTER OF MASS X-VEL     CENTER OF MASS Y-VEL     CENTER OF MASS Z-VEL      MAXIMUM TEMPERATURE          MAXIMUM DENSITY        MAXIMUM T_S / T_E
            0       0.0000000000000000   1.4902024346464787e+20   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   2.9583337617665660e+37   2.9583337617665660e+37   0.0000000000000000e+00   2.9583337617665660e+37   2.7306640614148750e+04   6.8067773166795473e+02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000156899e+08   3.1017131702472486e+07   0.0000000000000000e+00
-           1   5.9694928185132600e-08   1.4902038799260451e+20  -9.8852813993167483e+18   7.0758681407801181e+23  -3.7077151117729030e+15  -7.8559831399077499e+18   7.9789428545426686e+19   1.9321854546506737e+28   2.7622571688141327e+29   2.9583392229687520e+37   2.9583392505913144e+37   0.0000000000000000e+00   2.9583392505913144e+37   2.7306640614147251e+04   6.8067725884159722e+02   0.0000000000000000e+00  -6.6335093690719216e-02   4.7482550784469004e+03  -2.4880589573803199e-05   1.0000972452748042e+08   3.1017431302318867e+07   1.9852339015703824e-12
-           2   1.2535934918877847e-07   1.4902056479398973e+20  -2.0759055840781603e+19   1.4862415462814482e+24  -1.6351019859135748e+16  -3.4644556495600996e+19   3.5187152213757382e+20   4.0584338401784393e+28   5.5823502529134965e+29   2.9583452181108320e+37   2.9583452739343134e+37   0.0000000000000000e+00   2.9583452739343134e+37   2.7306640614140473e+04   6.8067697458181988e+02   0.0000000000000000e+00  -1.3930329595435043e-01   9.9733989623249036e+03  -1.0972324445112568e-04   1.0003863648020169e+08   3.1017667549506564e+07   1.9917944188178796e-12
+           1   5.9694928185132600e-08   1.4902038799260451e+20  -9.8852814000373391e+18   7.0758681407801342e+23  -3.7077151119947915e+15  -7.8559831403383900e+18   7.9789428550115738e+19   1.9321854546506743e+28   2.7622571688115688e+29   2.9583392229687520e+37   2.9583392505913144e+37   0.0000000000000000e+00   2.9583392505913144e+37   2.7306640614147251e+04   6.8067725884159722e+02   0.0000000000000000e+00  -6.6335093695554737e-02   4.7482550784469113e+03  -2.4880589575292177e-05   1.0000972452758998e+08   3.1017431302318867e+07   2.1834522208251888e-12
+           2   1.2535934918877847e-07   1.4902056479398969e+20  -2.0759055857352421e+19   1.4862415462814437e+24  -1.6351019861504854e+16  -3.4644556500407476e+19   3.5187152219157529e+20   4.0584338401784253e+28   5.5823502529036681e+29   2.9583452181108315e+37   2.9583452739343125e+37   0.0000000000000000e+00   2.9583452739343125e+37   2.7306640614140484e+04   6.8067697458182033e+02   0.0000000000000000e+00  -1.3930329606554864e-01   9.9733989623248763e+03  -1.0972324446702355e-04   1.0003863648118678e+08   3.1017667549506564e+07   1.8217100523183610e-12
diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out
@@ -2,4 +2,4 @@
 #   TIMESTEP                     TIME                 Mass He4                 Mass C12                 Mass O16                Mass Ne20                Mass Mg24                Mass Si28                 Mass S32                Mass Ar36                Mass Ca40                Mass Ti44                Mass Cr48                Mass Fe52                Mass Ni56
            0       0.0000000000000000   3.5899765861250859e-15   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.4944801491560746e-24   7.1354824912929246e-14
            1   5.9694928185132600e-08   3.5899765861251656e-15   7.4944874211898590e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.4944874177116286e-24   7.1354897598484802e-14
-           2   1.2535934918877847e-07   3.5899765861253060e-15   7.4944963166333664e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.1354986514892610e-14
+           2   1.2535934918877847e-07   3.5899765861253060e-15   7.4944963166333664e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.4944963093524191e-24   7.1354986514892597e-14
diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out
@@ -1,29 +1,29 @@
  plotfile = plt00095
  time = 1.25
               variables          minimum value          maximum value
- density                      8.7136176158e-05           13348283.786
- xmom                        -4.4636648292e+14       1.4969028735e+15
- ymom                        -1.8931343553e+15       1.9807937913e+15
+ density                      8.7136170328e-05           13347626.831
+ xmom                        -4.4623942798e+14       1.4982725823e+15
+ ymom                        -1.8879302767e+15        1.980304858e+15
  zmom                                        0                      0
- rho_E                        7.7947390767e+11       5.6568077669e+24
- rho_e                        7.4321344551e+11       5.6343409475e+24
- Temp                                   100000           3972527783.9
- rho_He4                      8.7136176158e-17           1473.9666386
- rho_C12                      3.4854470463e-05           5030539.1488
- rho_O16                      5.2281705694e-05           7778301.6377
- rho_Ne20                     8.7136176158e-17           1023673.1153
- rho_Mg24                     8.7136176158e-17           1040419.2782
- rho_Si28                     8.7136176158e-17           4251082.0739
- rho_S32                      8.7136176158e-17           2179431.2961
- rho_Ar36                     8.7136176158e-17           497747.48798
- rho_Ca40                     8.7136176158e-17             382056.037
- rho_Ti44                     8.7136176158e-17           1576.0930955
- rho_Cr48                     8.7136176158e-17           1467.9139369
- rho_Fe52                     8.7136176158e-17           14831.710059
- rho_Ni56                     8.7136176158e-17           182702.27304
- phiGrav                     -4.6147467267e+17      -2.2055818332e+16
- grav_x                          -461195258.85          -48603.568291
- grav_y                          -444709392.81           392306861.64
+ rho_E                        7.7948360809e+11         5.64675518e+24
+ rho_e                        7.4322366401e+11       5.6244082614e+24
+ Temp                                   100000           3969215375.1
+ rho_He4                      8.7136170328e-17           1473.6068478
+ rho_C12                      3.4854468131e-05           5028784.2519
+ rho_O16                      5.2281702196e-05           7774667.7686
+ rho_Ne20                     8.7136170328e-17           448739.14288
+ rho_Mg24                     8.7136170328e-17           1131012.4412
+ rho_Si28                     8.7136170328e-17           4244213.8084
+ rho_S32                      8.7136170328e-17           2178795.6667
+ rho_Ar36                     8.7136170328e-17           496791.83986
+ rho_Ca40                     8.7136170328e-17           381536.89513
+ rho_Ti44                     8.7136170328e-17           1576.1652647
+ rho_Cr48                     8.7136170328e-17           1467.4931808
+ rho_Fe52                     8.7136170328e-17           14841.830382
+ rho_Ni56                     8.7136170328e-17           182905.91385
+ phiGrav                     -4.6146743474e+17      -2.2055817929e+16
+ grav_x                          -461199998.92           -48603.52543
+ grav_y                          -444710966.19           392312565.37
  grav_z                                      0                      0
- rho_enuc                    -7.6356851771e+21       5.7259582003e+26
+ rho_enuc                    -7.4954583352e+21       6.7150736332e+26
 
diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H
@@ -64,6 +64,13 @@ 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)
+/// if r < model::profile(model_index).r(0) then we return 0, likewise
+/// if r > model::profile(model_index).r(model::npts-2) then we return model::npt-2,
+/// since this will give us the interval [npts-2, npts-1] to interpolate in
+///
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
 int
 locate(const Real r, const int model_index) {
@@ -74,7 +81,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 {
 
@@ -91,7 +98,7 @@ locate(const Real r, const int model_index) {
             }
         }
 
-        loc = ihi;
+        loc = ilo;
     }
 
     return loc;
@@ -102,73 +109,38 @@ 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.
+    // this gives us an index such that profile.r(id) < r < profile.r(id+1)
 
     int id = locate(r, model_index);
 
     Real slope;
     Real interp;
 
-    if (id == 0) {
-
-       slope = (model::profile(model_index).state(id+1, var_index) -
-                model::profile(model_index).state(id, var_index)) /
-           (model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
-       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),
-                                model::profile(model_index).state(id, var_index));
-       Real maxvar = amrex::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);
-
-    } else if (id == model::npts-1) {
-
-       slope = (model::profile(model_index).state(id, var_index) -
-                model::profile(model_index).state(id-1, var_index)) /
-           (model::profile(model_index).r(id) - model::profile(model_index).r(id-1));
-       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),
-                                model::profile(model_index).state(id, var_index));
-       Real maxvar = amrex::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);
-
-    } else {
-
-        if (r >= model::profile(model_index).r(id)) {
-
-            slope = (model::profile(model_index).state(id+1, var_index) -
-                     model::profile(model_index).state(id, var_index)) /
-                (model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
-            interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);
-
-        } else {
-
-            slope = (model::profile(model_index).state(id, var_index) -
-                     model::profile(model_index).state(id-1, var_index)) /
-                (model::profile(model_index).r(id) - model::profile(model_index).r(id-1));
-            interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);
-
-        }
-
+    slope = (model::profile(model_index).state(id+1, var_index) -
+             model::profile(model_index).state(id, var_index)) /
+        (model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
+    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.  We don't
+    // do this at the lower boundary, which usually corresponds to the center of the star.
+    if (r >= model::profile(model_index).r(0)) {
+	Real minvar = std::min(model::profile(model_index).state(id+1, var_index),
+			       model::profile(model_index).state(id, var_index));
+	Real maxvar = std::max(model::profile(model_index).state(id+1, var_index),
+			       model::profile(model_index).state(id, var_index));
+	interp = std::clamp(interp, minvar, maxvar);
     }
 
     return interp;
 
 }
 
-// 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)
 {

diff --git a/Util/model_parser/test/GNUmakefile b/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
diff --git a/Util/model_parser/test/Make.package b/Util/model_parser/test/Make.package
@@ -0,0 +1,4 @@
+CEXE_sources += main.cpp
+CEXE_sources += extern_parameters.cpp
+CEXE_headers += extern_parameters.H
+
diff --git a/Util/model_parser/test/README.md b/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.
diff --git a/Util/model_parser/test/main.cpp b/Util/model_parser/test/main.cpp
@@ -0,0 +1,69 @@
+#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);
+
+    // test locate and ensure that the index we get is such that
+    // profile.r(idx) < r < profile.r(idx+1)
+
+    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
+    // test to make sure that we see that
+
+    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));
+
+    // test the bounds.  Our model spans r = [5.e5, 4.0955e9]
+
+    std::cout << "testing boundaries of the model" << std::endl;
+
+    auto idx_lo_bnd = locate(-1.0, 0);
+    AMREX_ALWAYS_ASSERT(idx_lo_bnd == 0);
+    std::cout << "value a r < 0 = " << interpolate(-1.0, model::idens) << std::endl;
+
+    auto idx_hi_bnd = locate(4.1e9, 0);
+    AMREX_ALWAYS_ASSERT(idx_hi_bnd == model::npts-2);
+    std::cout << "value a r > r_max = " << interpolate(4.1e9, model::idens) << std::endl;
+
+    // test if we interpolate to exactly a point in the model that we
+    // recover the data there to roundoff
+
+    std::cout << "testing interpolating exactly on model point" << std::endl;
+
+    int idx_test{100};
+    Real r_test = model::profile(0).r(idx_test);
+    auto dens_test = interpolate(r_test, model::idens);
+
+    AMREX_ALWAYS_ASSERT(std::abs(dens_test - model::profile(0).state(idx_test, model::idens)) < 1.e-15_rt);
+
+}