Merge branch 'fix_interp' into add_cubic_interp

.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

Exec/science/massive_star/_prob_params

@@ -1,7 +1,10 @@
 model_name      character       ""        y
 
+interpolate_pres   int          0         y
+
 perturb_model   int             0         y
 velpert_amplitude  real         1.e5      y
 velpert_scale   real            1.e7      y
 
 isi28           int            -1
+

Exec/science/massive_star/problem_initialize_state_data.H

@@ -33,6 +33,7 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,URHO) = interpolate(dist, model::idens);
     state(i,j,k,UTEMP) = interpolate(dist, model::itemp);
+    Real pres = interpolate(dist, model::ipres);
     for (int n = 0; n < NumSpec; n++) {
         state(i,j,k,UFS+n) = amrex::max(interpolate(dist, model::ispec+n), small_x);
     }
@@ -116,22 +117,29 @@ void problem_initialize_state_data (int i, int j, int k,
         }
     }
 
-    // now call the EOS -- reload the eos_state since composition may
-    // have changed
+    // now call the EOS -- we need to use an eos_t here since we are
+    // going to take pressure as an input
 
-    burn_state.rho = state(i,j,k,URHO);
-    burn_state.T = state(i,j,k,UTEMP);
+    eos_t eos_state;
+    eos_state.rho = state(i,j,k,URHO);
+    eos_state.T = state(i,j,k,UTEMP);
+    eos_state.p = pres;
     for (int n = 0; n < NumSpec; n++) {
-        burn_state.xn[n] = state(i,j,k,UFS+n);
+        eos_state.xn[n] = state(i,j,k,UFS+n);
     }
     for (int n = 0; n < NumAux; n++) {
-        burn_state.aux[n] = state(i,j,k,UFX+n);
+        eos_state.aux[n] = state(i,j,k,UFX+n);
     }
 
-    eos(eos_input_rt, burn_state);
+    if (problem::interpolate_pres == 1) {
+        eos(eos_input_rp, eos_state);
+        state(i,j,k,UTEMP) = eos_state.T;
+    } else {
+        eos(eos_input_rt, eos_state);
+    }
 
-    state(i,j,k,UEINT) = state(i,j,k,URHO) * burn_state.e;
-    state(i,j,k,UEDEN) = state(i,j,k,URHO) * burn_state.e;
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,URHO) * eos_state.e;
 
     // add density scaling to composition
 

Source/driver/_cpp_parameters

@@ -32,20 +32,6 @@ update_sources_after_reflux  int           1
 allow_non_unit_aspect_zones  int           0
 
 
-#-----------------------------------------------------------------------------
-# category: initialization
-#-----------------------------------------------------------------------------
-
-# for fourth order, we usually assume that the initialization is done
-# to cell centers and we convert to cell-averages.  With this option,
-# we take the initialization as cell-averages (except for T, which we
-# compute to fourth-order through the EOS after initialization).
-initialization_is_cell_average  int        0
-
-# type of interpolation to use for mapping the 1D initial model onto the
-# domain.  0 = linear, 1 = cubic
-model_interpolation_method      int        0
-
 #-----------------------------------------------------------------------------
 # category: hydrodynamics
 #-----------------------------------------------------------------------------
@@ -76,6 +62,12 @@ time_integration_method      int           0
 # do we use a limiter with the fourth-order accurate reconstruction?
 limit_fourth_order           int           1
 
+# for fourth order, we usually assume that the initialization is done
+# to cell centers and we convert to cell-averages.  With this option,
+# we take the initialization as cell-averages (except for T, which we
+# compute to fourth-order through the EOS after initialization).
+initialization_is_cell_average  int        0
+
 # should we use a reconstructed version of Gamma_1 in the Riemann
 # solver? or the default zone average (requires SDC
 # integration, since we do not trace)

Util/model_parser/model_parser.H

@@ -65,8 +65,8 @@ namespace model_string
 
 
 ///
-/// return the index into the model coordinate such that
-/// model::profile(model_index).r(lo) < r < model::profile(model_index).r(hi)
+/// 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
@@ -98,13 +98,13 @@ locate(const Real r, const int model_index) {
         loc = ihi;
     }
 
-    return loc;
+    return loc-1;
 }
 
 
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
 Real
-linear_interpolate(const Real r, const int var_index, const int model_index=0) {
+interpolate(const Real r, const int var_index, const int model_index=0) {
 
     int id = locate(r, model_index);
 
@@ -165,68 +165,6 @@ linear_interpolate(const Real r, const int var_index, const int model_index=0) {
 }
 
 
-AMREX_INLINE AMREX_GPU_HOST_DEVICE
-Real
-cubic_interpolate(const Real r, const int var_index, const int model_index=0) {
-
-    int id = locate(r, model_index);
-    // we will use 4 zones, id-1, id, id+1, id+2
-    // make sure all are valid indices
-    id = std::clamp(id, 1, model::npts-3);
-
-    // note: we are assuming that everything is equally spaced here
-
-    // fit a cubic of the form
-    // v(r) = a (r - r_i)**3 + b (r - r_i)**2 + c (r - r_i) + d
-    // to the data (rs, vs)
-    // we take r_i to be r[1]
-
-    const Real vs[4] = {model::profile(model_index).state(id-1, var_index),
-                        model::profile(model_index).state(id, var_index),
-                        model::profile(model_index).state(id+1, var_index),
-                        model::profile(model_index).state(id+2, var_index)};
-
-    const Real rs[4] = {model::profile(model_index).r(id-1),
-                        model::profile(model_index).r(id),
-                        model::profile(model_index).r(id+1),
-                        model::profile(model_index).r(id+2)};
-
-    const Real dx = rs[1] - rs[0];
-
-    Real a = (3 * vs[1] - 3 * vs[2] + vs[3] - vs[0]) / (6 * std::pow(dx, 3));
-    Real b = (-2 * vs[1] + vs[2] + vs[0]) / (2 * dx * dx);
-    Real c = (-3 * vs[1] + 6 * vs[2] - vs[3] - 2 * vs[0]) / (6 * dx);
-    Real d = vs[1];
-
-    Real interp = a * std::pow(r - rs[1], 3) + b * std::pow(r - rs[1], 2) + c * (r - rs[1]) + d;
-
-    // safety check to make sure interp lies within the bounding points
-    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;
-}
-
-
-///
-/// Find the value of model_state component var_index at point r.
-/// This is a wrapper for the specific implementation of the interpolation.
-///
-AMREX_INLINE AMREX_GPU_HOST_DEVICE
-Real
-interpolate(const Real r, const int var_index, const int model_index=0) {
-
-    if (model_interpolation_method == 0) {
-        return linear_interpolate(r, var_index, model_index);
-    } else {
-        return cubic_interpolate(r, var_index, model_index);
-    }
-}
-
-
 ///
 /// 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.

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));
+
+}