don't limit for r < r(0)

AMReX-Astro · Dec 22, 2023 · 8cc53b0 · 8cc53b0
1 parent db0fedf
commit 8cc53b0
Showing 1 changed file with 9 additions and 6 deletions.
diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H
@@ -122,12 +122,15 @@ 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 = 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);
+    // 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;