merge

Source/driver/_cpp_parameters

@@ -32,6 +32,20 @@ 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
 #-----------------------------------------------------------------------------
@@ -62,12 +76,6 @@ 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,7 +65,7 @@ namespace model_string
 
 
 ///
-/// return the index into the model coordinate, loc, such that
+/// return the index into the model coordinate such that
 /// model::profile(model_index).r(loc) < r < model::profile(model_index).r(loc+1)
 ///
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
@@ -104,7 +104,7 @@ locate(const Real r, const int model_index) {
 
 AMREX_INLINE AMREX_GPU_HOST_DEVICE
 Real
-interpolate(const Real r, const int var_index, const int model_index=0) {
+linear_interpolate(const Real r, const int var_index, const int model_index=0) {
 
     int id = locate(r, model_index);
 
@@ -165,6 +165,68 @@ 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.