diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 774d6916ac..5374500539 100644 --- a/Source/driver/_cpp_parameters +++ b/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) diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index e8cd0c6f5f..64b7035585 100644 --- a/Util/model_parser/model_parser.H +++ b/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.