From 362b7b754c3e71b767d5d1ba1357ba3b3ddbb377 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 30 Nov 2023 12:31:46 -0500 Subject: [PATCH] fix an offset in the interpolation --- .../science/massive_star/problem_initialize.H | 1 - Source/driver/_cpp_parameters | 20 ++---- Util/model_parser/model_parser.H | 70 ++----------------- 3 files changed, 10 insertions(+), 81 deletions(-) diff --git a/Exec/science/massive_star/problem_initialize.H b/Exec/science/massive_star/problem_initialize.H index 6e44398280..c2fd00cf88 100644 --- a/Exec/science/massive_star/problem_initialize.H +++ b/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; diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 5374500539..774d6916ac 100644 --- a/Source/driver/_cpp_parameters +++ b/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) diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index b4d25dcf62..e8cd0c6f5f 100644 --- a/Util/model_parser/model_parser.H +++ b/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.