From 1f01da8e62850fdfc0d9bcd73b0467a8824d7c8d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 4 Jan 2024 10:17:34 -0500 Subject: [PATCH 1/3] constexpr the runtime parameters --- Exec/Make.auto_source | 2 +- Source/driver/Castro.cpp | 15 +++-- Source/driver/Castro_advance_ctu.cpp | 4 +- Source/driver/Castro_io.cpp | 6 +- Source/driver/Castro_setup.cpp | 14 ++-- Source/driver/parse_castro_params.py | 99 ++++++++++++++++------------ Source/gravity/Castro_pointmass.cpp | 2 +- Source/gravity/Gravity.cpp | 2 +- Source/scf/scf_relax.cpp | 8 +-- 9 files changed, 88 insertions(+), 64 deletions(-) diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source index 178e0f43e2..77b48e8fb3 100644 --- a/Exec/Make.auto_source +++ b/Exec/Make.auto_source @@ -69,7 +69,7 @@ CPP_PARAMETERS := $(TOP)/Source/driver/_cpp_parameters $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H: $(CPP_PARAMETERS) @if [ ! -d $(CASTRO_AUTO_SOURCE_DIR) ]; then mkdir -p $(CASTRO_AUTO_SOURCE_DIR); fi - PYTHONPATH=$(MICROPHYSICS_HOME)/util/build_scripts $(TOP)/Source/driver/parse_castro_params.py -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS) + PYTHONPATH=$(MICROPHYSICS_HOME)/util/build_scripts $(TOP)/Source/driver/parse_castro_params.py --constexpr -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS) # for debugging test_cxx_params: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index f10c59b0e7..ae47a6b380 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -415,9 +415,9 @@ Castro::read_params () } // Make sure not to call refluxing if we're not actually doing any hydro. - if (do_hydro == 0) { - do_reflux = 0; - } +// if (do_hydro == 0) { +// do_reflux = 0; +// } if (max_dt < fixed_dt) { @@ -475,9 +475,9 @@ Castro::read_params () } } - if (dgeom.IsRZ() == 1) { - rot_axis = 2; - } +// if (dgeom.IsRZ() == 1) { +// rot_axis = 2; +// } #if (AMREX_SPACEDIM == 1) if (do_rotation == 1) { std::cerr << "ERROR:Castro::Rotation not implemented in 1d\n"; @@ -2070,7 +2070,8 @@ Castro::post_timestep (int iteration_local) Real max_field_val = std::numeric_limits::min(); for (int lev = 0; lev <= parent->finestLevel(); ++lev) { - auto mf = getLevel(lev).derive(castro::stopping_criterion_field, state[State_Type].curTime(), 0); + std::string scf{castro::stopping_criterion_field}; + auto mf = getLevel(lev).derive(scf, state[State_Type].curTime(), 0); max_field_val = std::max(max_field_val, mf->max(0)); } diff --git a/Source/driver/Castro_advance_ctu.cpp b/Source/driver/Castro_advance_ctu.cpp index 911a2780e9..ebd22a4480 100644 --- a/Source/driver/Castro_advance_ctu.cpp +++ b/Source/driver/Castro_advance_ctu.cpp @@ -364,7 +364,7 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // approximately compensate for that. if (in_retry) { - sdc_iters += 1; + //sdc_iters += 1; amrex::Print() << "Adding an SDC iteration due to the retry." << std::endl << std::endl; } @@ -412,7 +412,7 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // Set sdc_iters to its original value, in case we modified it above. - sdc_iters = sdc_iters_old; + //sdc_iters = sdc_iters_old; // If we're allowing for retries, check for that here. diff --git a/Source/driver/Castro_io.cpp b/Source/driver/Castro_io.cpp index ea0336669e..e31697d60c 100644 --- a/Source/driver/Castro_io.cpp +++ b/Source/driver/Castro_io.cpp @@ -181,7 +181,8 @@ Castro::restart (Amr& papa, PMFile.open(FullPathPMFile.c_str(), std::ios::in); if (PMFile.good()) { - PMFile >> point_mass; + Real pm; + PMFile >> pm; // point_mass; PMFile.close(); } @@ -198,7 +199,8 @@ Castro::restart (Amr& papa, RotationFile.open(FullPathRotationFile.c_str(), std::ios::in); if (RotationFile.is_open()) { - RotationFile >> castro::rotational_period; + Real rp; + RotationFile >> rp; //castro::rotational_period; amrex::Print() << " Based on the checkpoint, setting the rotational period to " << std::setprecision(7) << std::fixed << castro::rotational_period << " s.\n"; RotationFile.close(); diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp index 9506fdde2c..4c22e10f03 100644 --- a/Source/driver/Castro_setup.cpp +++ b/Source/driver/Castro_setup.cpp @@ -209,6 +209,7 @@ Castro::variableSetUp () // set small positive values of the "small" quantities if they are // negative +#if 0 if (small_dens < 0.0_rt) { small_dens = 1.e-100_rt; } @@ -224,13 +225,16 @@ Castro::variableSetUp () if (small_ener < 0.0_rt) { small_ener = 1.e-100_rt; } +#endif // now initialize the C++ Microphysics #ifdef REACTIONS network_init(); #endif - eos_init(castro::small_temp, castro::small_dens); + Real st{castro::small_temp}; + Real sd{castro::small_dens}; + eos_init(st, sd); #ifdef RADIATION opacity_init(); @@ -240,11 +244,11 @@ Castro::variableSetUp () // with the minimum permitted by the EOS, and vice versa. Real new_min_T = std::max(small_temp, EOSData::mintemp); - small_temp = new_min_T; + //small_temp = new_min_T; EOSData::mintemp = new_min_T; Real new_min_rho = std::max(small_dens, EOSData::mindens); - small_dens = new_min_rho; + //small_dens = new_min_rho; EOSData::mindens = new_min_rho; // Given small_temp and small_dens, compute small_pres @@ -267,8 +271,8 @@ Castro::variableSetUp () eos(eos_input_rt, eos_state); - castro::small_pres = amrex::max(castro::small_pres, eos_state.p); - castro::small_ener = amrex::max(castro::small_ener, eos_state.e); + //castro::small_pres = amrex::max(castro::small_pres, eos_state.p); + //castro::small_ener = amrex::max(castro::small_ener, eos_state.e); // some consistency checks on the parameters #ifdef REACTIONS diff --git a/Source/driver/parse_castro_params.py b/Source/driver/parse_castro_params.py index aedeea6624..f742be2bba 100755 --- a/Source/driver/parse_castro_params.py +++ b/Source/driver/parse_castro_params.py @@ -129,7 +129,7 @@ def read_param_file(infile): return params -def write_headers(params, out_directory, struct_name): +def write_headers(params, out_directory, struct_name, constexpr): # output @@ -152,19 +152,23 @@ def write_headers(params, out_directory, struct_name): cd.write(f"#ifndef {nm.upper()}_DECLARES_H\n") cd.write(f"#define {nm.upper()}_DECLARES_H\n") - cd.write("\n") - cd.write(f"namespace {nm} {{\n") - - for ifdef in ifdefs: - if ifdef is None: - for p in [q for q in params_nm if q.ifdef is None]: - cd.write(p.get_declare_string()) - else: - cd.write(f"#ifdef {ifdef}\n") - for p in [q for q in params_nm if q.ifdef == ifdef]: - cd.write(p.get_declare_string()) - cd.write("#endif\n") - cd.write("}\n\n") + if constexpr: + pass + else: + cd.write("\n") + cd.write(f"namespace {nm} {{\n") + + for ifdef in ifdefs: + if ifdef is None: + for p in [q for q in params_nm if q.ifdef is None]: + cd.write(p.get_declare_string()) + else: + cd.write(f"#ifdef {ifdef}\n") + for p in [q for q in params_nm if q.ifdef == ifdef]: + cd.write(p.get_declare_string()) + cd.write("#endif\n") + + cd.write("}\n\n") cd.write("#endif\n") cd.close() @@ -177,19 +181,27 @@ def write_headers(params, out_directory, struct_name): cp.write(CWARNING) cp.write(f"#ifndef {nm.upper()}_PARAMS_H\n") cp.write(f"#define {nm.upper()}_PARAMS_H\n") - + cp.write("#include \n") + cp.write(f"using namespace amrex::literals;\n") cp.write("\n") cp.write(f"namespace {nm} {{\n") - for ifdef in ifdefs: - if ifdef is None: - for p in [q for q in params_nm if q.ifdef is None]: - cp.write(p.get_declare_string(with_extern=True)) - else: - cp.write(f"#ifdef {ifdef}\n") - for p in [q for q in params_nm if q.ifdef == ifdef]: - cp.write(p.get_declare_string(with_extern=True)) - cp.write("#endif\n") + if constexpr: + for p in params_nm: + type = p.get_cxx_decl() + if type == "std::string": + type = "std::string_view" + cp.write(f"constexpr {type} {p.cpp_var_name}{{{p.default_format()}}};\n") + else: + for ifdef in ifdefs: + if ifdef is None: + for p in [q for q in params_nm if q.ifdef is None]: + cp.write(p.get_declare_string(with_extern=True)) + else: + cp.write(f"#ifdef {ifdef}\n") + for p in [q for q in params_nm if q.ifdef == ifdef]: + cp.write(p.get_declare_string(with_extern=True)) + cp.write("#endif\n") cp.write("}\n\n") cp.write("#endif\n") cp.close() @@ -202,22 +214,25 @@ def write_headers(params, out_directory, struct_name): cq.write(CWARNING) - for ifdef in ifdefs: - if ifdef is None: - for p in [q for q in params_nm if q.ifdef is None]: - cq.write(p.get_default_string()) - cq.write(p.get_query_string()) - cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro")) - cq.write("\n") - else: - cq.write(f"#ifdef {ifdef}\n") - for p in [q for q in params_nm if q.ifdef == ifdef]: - cq.write(p.get_default_string()) - cq.write(p.get_query_string()) - cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro")) - cq.write("\n") - cq.write("#endif\n") - cq.write("\n") + if constexpr: + pass + else: + for ifdef in ifdefs: + if ifdef is None: + for p in [q for q in params_nm if q.ifdef is None]: + cq.write(p.get_default_string()) + cq.write(p.get_query_string()) + cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro")) + cq.write("\n") + else: + cq.write(f"#ifdef {ifdef}\n") + for p in [q for q in params_nm if q.ifdef == ifdef]: + cq.write(p.get_default_string()) + cq.write(p.get_query_string()) + cq.write(p.get_query_struct_string(struct_name=struct_name, class_name="Castro")) + cq.write("\n") + cq.write("#endif\n") + cq.write("\n") cq.close() # write the job info tests @@ -289,13 +304,15 @@ def main(): help="output directory for the generated files") parser.add_argument("-s", type=str, default="params", help="name for the name struct that will hold the parameters") + parser.add_argument("--constexpr", action="store_true", + help="force the parameters to be constexpr without being able to query") parser.add_argument("input_file", type=str, nargs=1, help="input file containing the list of parameters we will define") args = parser.parse_args() p = read_param_file(args.input_file[0]) - write_headers(p, args.o, args.s) + write_headers(p, args.o, args.s, args.constexpr) if __name__ == "__main__": main() diff --git a/Source/gravity/Castro_pointmass.cpp b/Source/gravity/Castro_pointmass.cpp index 1bbdca2050..485a3b6886 100644 --- a/Source/gravity/Castro_pointmass.cpp +++ b/Source/gravity/Castro_pointmass.cpp @@ -103,7 +103,7 @@ Castro::pointmass_update(Real time, Real dt) << ", to " << point_mass + mass_change_at_center << std::endl << std::endl; } - point_mass += mass_change_at_center; + //point_mass += mass_change_at_center; #ifdef _OPENMP #pragma omp parallel diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 6665b0ef38..df6b046211 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -342,7 +342,7 @@ Gravity::install_level (int level, std::string Gravity::get_gravity_type() { - return gravity::gravity_type; + return std::string{gravity::gravity_type}; } int Gravity::get_max_solve_level() diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp index 1add2074e3..796a188a0b 100644 --- a/Source/scf/scf_relax.cpp +++ b/Source/scf/scf_relax.cpp @@ -327,7 +327,7 @@ Castro::do_hscf_solve() if (std::abs(scf_polar_radius - scf_equatorial_radius) / std::abs(scf_equatorial_radius) < 1.e-6) { - rotational_period = 0.0; + //rotational_period = 0.0; } else { @@ -344,9 +344,9 @@ Castro::do_hscf_solve() // Let's also be sure not to let the period // change by too much in a single iteration. - rotational_period = amrex::min(1.1_rt * rotational_period, - amrex::max(0.9_rt * rotational_period, - 2.0_rt * M_PI / omega)); + //rotational_period = amrex::min(1.1_rt * rotational_period, + // amrex::max(0.9_rt * rotational_period, + // 2.0_rt * M_PI / omega)); } From 03761447f6378e48d17c863825e815caa5b2f53e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 4 Jan 2024 11:05:13 -0500 Subject: [PATCH 2/3] hardcode the parameters we use for subchandra --- Source/driver/_cpp_parameters | 54 +++++++++++++++++------------------ 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 4d1dde2581..90cec3c607 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -25,7 +25,7 @@ do_reflux int 1 # whether to re-compute new-time source terms after a reflux # Note: this only works for the CTU and simple-SDC time_integration_method # drivers -update_sources_after_reflux int 1 +update_sources_after_reflux int 0 # Castro was originally written assuming dx = dy = dz. This assumption is # enforced at runtime. Setting allow_non_unit_aspect_zones = 1 opts out. @@ -40,24 +40,24 @@ allow_non_unit_aspect_zones int 0 difmag Real 0.1 # the small density cutoff. Densities below this value will be reset -small_dens Real -1.e200 +small_dens Real 1.e-5 # the small temperature cutoff. Temperatures below this value will # be reset -small_temp Real -1.e200 +small_temp Real 1.e5 # the small pressure cutoff. Pressures below this value will be reset -small_pres Real -1.e200 +small_pres Real 5.0585e7 # the small specific internal energy cutoff. Internal energies below this # value will be reset -small_ener Real -1.e200 +small_ener Real 7.49081e+12 # permits hydro to be turned on and off for running pure rad problems -do_hydro int -1 +do_hydro int 1 # how do we advance in time? 0 = CTU + Strang, 1 is not used, 2 = SDC, 3 = simplified-SDC -time_integration_method int 0 +time_integration_method int 3 # do we use a limiter with the fourth-order accurate reconstruction? limit_fourth_order int 1 @@ -112,7 +112,7 @@ hybrid_riemann int 0 # 0: Colella, Glaz, \& Ferguson (a two-shock solver); # 1: Colella \& Glaz (a two-shock solver) # 2: HLLC -riemann_solver int 0 +riemann_solver int 1 # for the Colella \& Glaz Riemann solver, the maximum number # of iterations to take when solving for the star state @@ -174,7 +174,7 @@ limit_fluxes_on_small_dens int 0 speed_limit Real 0.0 # permits sponge to be turned on and off -do_sponge int 0 +do_sponge int 1 # if we are using the sponge, whether to use the implicit solve for it sponge_implicit int 1 @@ -312,16 +312,16 @@ max_dt Real 1.e200 # the effective Courant number to use---we will not allow the hydrodynamic # waves to cross more than this fraction of a zone over a single timestep -cfl Real 0.8 +cfl Real 0.2 # a factor by which to reduce the first timestep from that requested by # the timestep estimators -init_shrink Real 1.0 +init_shrink Real 0.05 # the maximum factor by which the timestep can increase or decrease from # one step to the next. Must be greater than 1.0---use max_dt to set a cap # on the timestep. -change_max Real 1.1 +change_max Real 1.025 # whether to check that we will take a valid timestep before the advance check_dt_before_advance int 1 @@ -355,14 +355,14 @@ abundance_failure_tolerance Real 1.e-2 # Do not abort for invalid species abundances if the zone's density is below # this threshold. -abundance_failure_rho_cutoff Real -1.e200 +abundance_failure_rho_cutoff Real 1 # Regrid after every timestep. use_post_step_regrid int 0 # Do not permit more subcycled timesteps than this parameter. # Set to a negative value to disable this criterion. -max_subcycles int 10 +max_subcycles int 32 # Number of iterations for the simplified SDC advance. sdc_iters int 2 @@ -393,16 +393,16 @@ dtnuc_X Real 1.e200 dtnuc_X_threshold Real 1.e-3 # permits reactions to be turned on and off -- mostly for efficiency's sake -do_react int -1 +do_react int 1 # minimum temperature for allowing reactions to occur in a zone -react_T_min Real 0.0 +react_T_min Real 5.e7 # maximum temperature for allowing reactions to occur in a zone react_T_max Real 1.e200 # minimum density for allowing reactions to occur in a zone -react_rho_min Real 0.0 +react_rho_min Real 1 # maximum density for allowing reactions to occur in a zone react_rho_max Real 1.e200 @@ -450,7 +450,7 @@ diffuse_cond_scale_fac Real 1.0 DIFFUSION #----------------------------------------------------------------------------- # permits gravity calculation to be turned on and off -do_grav int -1 +do_grav int 1 # to we recompute the center used for the multipole gravity solve each step? moving_center int 0 @@ -518,10 +518,10 @@ sponge_lower_radius Real -1.0 SPONGE sponge_upper_radius Real -1.0 SPONGE # Minimum density at which to start applying the sponge -sponge_lower_density Real -1.0 SPONGE +sponge_lower_density Real 100 SPONGE # Density at which the sponge is fully applied -sponge_upper_density Real -1.0 SPONGE +sponge_upper_density Real 10000 SPONGE # Minimum pressure at which to start applying the sponge sponge_lower_pressure Real -1.0 SPONGE @@ -545,7 +545,7 @@ sponge_target_y_velocity Real 0.0 SPONGE sponge_target_z_velocity Real 0.0 SPONGE # Timescale on which the sponge operates -sponge_timescale Real -1.0 SPONGE +sponge_timescale Real 0.001 SPONGE #----------------------------------------------------------------------------- @@ -606,7 +606,7 @@ max_tagging_radius real 10.0e0 #----------------------------------------------------------------------------- # verbosity level (higher numbers mean more output) -(v, verbose) int 0 +(v, verbose) int 1 # do we dump the old state into the checkpoint files too? dump_old bool false @@ -620,7 +620,7 @@ domain_is_plane_parallel bool false print_update_diagnostics int (0, 1) # how often (number of coarse timesteps) to compute integral sums (for runtime diagnostics) -sum_interval int -1 +sum_interval int 5 # how often (simulation time) to compute integral sums (for runtime diagnostics) sum_per Real -1.0e0 @@ -650,7 +650,7 @@ store_omegadot int 0 # Do we store the burn weights as a diagnostic in the plotfile? Note, if this option is # enabled then more memory will be allocated to hold the results of the burn -store_burn_weights int 0 +store_burn_weights int 1 # Do we abort the run if the inputs file specifies a runtime parameter that we don't # know about? Note: this will only take effect for those namespaces where 100% @@ -705,7 +705,7 @@ timestamp_temperature int 0 @namespace: gravity # what type -gravity_type string "fillme" +gravity_type string "PoissonGrav" # if doing constant gravity, what is the acceleration @@ -720,7 +720,7 @@ drdxfac int 1 # the maximum mulitpole order to use for multipole BCs when doing # Poisson gravity -(max_multipole_order, lnum) int 0 +(max_multipole_order, lnum) int 6 # the level of verbosity for the gravity solve (higher number means more # output on the status of the solve / multigrid @@ -733,7 +733,7 @@ no_sync int 0 # gets us closer to the composite solution? This makes the # resulting fine grid calculation slightly more accurate, # at the cost of an additional Poisson solve per timestep. -do_composite_phi_correction int 1 +do_composite_phi_correction int 0 # For all gravity types, we can choose a maximum level for explicitly # calculating the gravity and associated potential. Above that level, From aa8277c3e54ff5245a179cfb093f087b0ef28859 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 4 Jan 2024 12:55:59 -0500 Subject: [PATCH 3/3] hook in Microphysics --- Exec/Make.auto_source | 2 +- Source/driver/Castro.cpp | 7 ------- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source index 77b48e8fb3..5c79e3ff5d 100644 --- a/Exec/Make.auto_source +++ b/Exec/Make.auto_source @@ -25,7 +25,7 @@ $(CASTRO_AUTO_SOURCE_DIR)/extern_parameters.H: $(EXTERN_PARAMETERS) @if [ ! -d $(CASTRO_AUTO_SOURCE_DIR) ]; then mkdir -p $(CASTRO_AUTO_SOURCE_DIR); fi $(MICROPHYSICS_HOME)/util/build_scripts/write_probin.py \ --cxx_prefix $(CASTRO_AUTO_SOURCE_DIR)/extern --use_namespace \ - --pa "$(EXTERN_PARAMETERS)" + --pa "$(EXTERN_PARAMETERS)" --constexpr #------------------------------------------------------------------------------ diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index ae47a6b380..02d683d09c 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -324,13 +324,6 @@ Castro::read_params () } #endif -#ifdef REACTIONS -#ifdef SIMPLIFIED_SDC - if (jacobian == 1) { - amrex::Abort("Simplified SDC requires the numerical Jacobian now (jacobian = 2)"); - } -#endif -#endif // sanity checks if (grown_factor < 1) {