fix an offset in the interpolation

Exec/science/massive_star/problem_initialize.H

@@ -19,7 +19,6 @@ void problem_initialize ()
 
     read_model_file(problem::model_name);
 
-
 #if AMREX_SPACEDIM == 1
     problem::center[0] = 0.0_rt;
 

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.