From fe81e5d363dea3c5eaaf05abe28f9b243b1c25c2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 26 Jun 2024 07:52:24 -0400 Subject: [PATCH 01/25] fix some namespace issues + add missing headers to runtime_parameters.H this will make it easier to use this for non-castro things (like exact_riemann) --- Source/driver/runtime_parameters.H | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/Source/driver/runtime_parameters.H b/Source/driver/runtime_parameters.H index 0619070811..4881668f2b 100644 --- a/Source/driver/runtime_parameters.H +++ b/Source/driver/runtime_parameters.H @@ -1,6 +1,7 @@ #ifndef RUNTIME_PARAMETERS_H #define RUNTIME_PARAMETERS_H +#include #include #include @@ -25,34 +26,34 @@ initialize_cpp_runparams() { { - ParmParse pp("castro"); + amrex::ParmParse pp("castro"); #include } #ifdef AMREX_PARTICLES { - ParmParse pp("particles"); + amrex::ParmParse pp("particles"); #include } #endif #ifdef DIFFUSION { - ParmParse pp("diffusion"); + amrex::ParmParse pp("diffusion"); #include } #endif #ifdef GRAVITY { - ParmParse pp("gravity"); + amrex::ParmParse pp("gravity"); #include } #endif #ifdef RADIATION { - ParmParse pp("radiation"); + amrex::ParmParse pp("radiation"); #include } #endif @@ -102,14 +103,14 @@ validate_runparams() for (const auto& nm: check_namespaces) { // "castro" - if (ParmParse::hasUnusedInputs(nm)) { + if (amrex::ParmParse::hasUnusedInputs(nm)) { amrex::Print() << "Warning: the following " + nm + ".* parameters are ignored\n"; - auto unused = ParmParse::getUnusedInputs(nm); + auto unused = amrex::ParmParse::getUnusedInputs(nm); for (const auto& p: unused) { amrex::Print() << p << "\n"; } amrex::Print() << std::endl; - if (abort_on_invalid_params) { + if (castro::abort_on_invalid_params) { amrex::Error("Error: invalid parameters"); } } From 000768fca4a2332245d320f32b90a57d4d882954 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 10:45:59 -0400 Subject: [PATCH 02/25] more work on getting castro params into exact Riemann --- Exec/Make.auto_source | 8 ++++- Source/driver/parse_castro_params.py | 13 +++++--- Util/exact_riemann/GNUmakefile | 7 ++++ Util/exact_riemann/main.cpp | 49 +++++++++++++++++++++------- 4 files changed, 60 insertions(+), 17 deletions(-) diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source index e9e93fc284..dc3dfe194b 100644 --- a/Exec/Make.auto_source +++ b/Exec/Make.auto_source @@ -72,9 +72,15 @@ CPP_PARAMETERS := $(TOP)/Source/driver/_cpp_parameters $(CASTRO_AUTO_SOURCE_DIR)/runtime_params.cpp: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H +STRUCT_USE_CASTRO_CLASS ?= TRUE +PARSE_ARGS := +ifneq ($(STRUCT_USE_CASTRO_CLASS), TRUE) + PARSE_ARGS += --without-castro-class +endif + $(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 $(PARSE_ARGS) -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS) # for debugging test_cxx_params: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H diff --git a/Source/driver/parse_castro_params.py b/Source/driver/parse_castro_params.py index 89f91a8ea7..4d869d09dd 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_and_source(params, out_directory, struct_name): +def write_headers_and_source(params, out_directory, struct_name, without_castro_class): # output @@ -176,19 +176,23 @@ def write_headers_and_source(params, out_directory, struct_name): cq.write(CWARNING) + class_name = "Castro" + if without_castro_class: + class_name = None + 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(p.get_query_struct_string(struct_name=struct_name, class_name=class_name)) 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(p.get_query_struct_string(struct_name=struct_name, class_name=class_name)) cq.write("\n") cq.write("#endif\n") cq.write("\n") @@ -305,13 +309,14 @@ 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("--without-castro-class", action="store_true", help="don't include Castro:: in the struct namespace") 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_and_source(p, args.o, args.s) + write_headers_and_source(p, args.o, args.s, args.without_castro_class) if __name__ == "__main__": main() diff --git a/Util/exact_riemann/GNUmakefile b/Util/exact_riemann/GNUmakefile index 76bc8e8a9f..8be2709e1a 100644 --- a/Util/exact_riemann/GNUmakefile +++ b/Util/exact_riemann/GNUmakefile @@ -23,9 +23,16 @@ EOS_DIR := helmholtz NETWORK_DIR := general_null NETWORK_INPUTS = ignition.net +# for the Castro runtime parameters, we don't want to use Castro:: +STRUCT_USE_CASTRO_CLASS := FALSE + EXTERN_SEARCH += . Bpack := ./Make.package Blocs := . +# we explicitly want runtime_parameters.H so we have access to +# Castro's runtime parameter system +Blocs += $(CASTRO_HOME)/Source/driver + include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp index c2b3ed296d..6e5a5f7150 100644 --- a/Util/exact_riemann/main.cpp +++ b/Util/exact_riemann/main.cpp @@ -13,27 +13,52 @@ #include -int main(int argc, char *argv[]) { +#include +#include +#include - amrex::Initialize(argc, argv); +// create a dummy Castro class just to hold the runtime parameters struct - std::cout << "starting the exact Riemann solver..." << std::endl; +class Castro { + + public: + + static params_t params; - // initialize the external runtime parameters in C++ + void doit() { - init_prob_parameters(); + initialize_cpp_runparams(); - init_extern_parameters(); + // initialize the external runtime parameters in C++ - // now initialize the C++ Microphysics + init_prob_parameters(); + + init_extern_parameters(); + + // now initialize the C++ Microphysics #ifdef REACTIONS - network_init(); + network_init(); #endif - amrex::Real small_temp = 1.e-200; - amrex::Real small_dens = 1.e-200; - eos_init(small_temp, small_dens); + amrex::Real small_temp = 1.e-200; + amrex::Real small_dens = 1.e-200; + eos_init(small_temp, small_dens); + + exact_riemann(); + + } + +}; + + +int main(int argc, char *argv[]) { + + amrex::Initialize(argc, argv); + + std::cout << "starting the exact Riemann solver..." << std::endl; - exact_riemann(); + // initialize the Castro runtime parameters + Castro castro; + castro.doit(); } From d01fa53e94e32ea453e73ef633bf7949088ad9f1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 12:32:33 -0400 Subject: [PATCH 03/25] this seems to work --- Util/exact_riemann/Make.package | 3 +++ Util/exact_riemann/main.cpp | 48 ++++++++++++--------------------- 2 files changed, 20 insertions(+), 31 deletions(-) diff --git a/Util/exact_riemann/Make.package b/Util/exact_riemann/Make.package index b823fe2404..2cc1befab7 100644 --- a/Util/exact_riemann/Make.package +++ b/Util/exact_riemann/Make.package @@ -6,3 +6,6 @@ CEXE_headers += riemann_sample.H CEXE_headers += riemann_support.H CEXE_sources += extern_parameters.cpp CEXE_headers += extern_parameters.H + +# automatically generated in tmp_build_dir +CEXE_sources += runtime_params.cpp \ No newline at end of file diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp index 6e5a5f7150..7f75f6fc36 100644 --- a/Util/exact_riemann/main.cpp +++ b/Util/exact_riemann/main.cpp @@ -14,51 +14,37 @@ #include #include -#include #include -// create a dummy Castro class just to hold the runtime parameters struct +// a global struct to hold the params +#include -class Castro { +int main(int argc, char *argv[]) { - public: + amrex::Initialize(argc, argv); - static params_t params; + std::cout << "starting the exact Riemann solver..." << std::endl; - void doit() { + // initialize the Castro runtime parameters - initialize_cpp_runparams(); + amrex::ParmParse pp("castro"); +#include - // initialize the external runtime parameters in C++ + // initialize the external runtime parameters in C++ - init_prob_parameters(); + init_prob_parameters(); - init_extern_parameters(); + init_extern_parameters(); - // now initialize the C++ Microphysics + // now initialize the C++ Microphysics #ifdef REACTIONS - network_init(); + network_init(); #endif - amrex::Real small_temp = 1.e-200; - amrex::Real small_dens = 1.e-200; - eos_init(small_temp, small_dens); - - exact_riemann(); - - } - -}; - + amrex::Real small_temp = 1.e-200; + amrex::Real small_dens = 1.e-200; + eos_init(small_temp, small_dens); -int main(int argc, char *argv[]) { - - amrex::Initialize(argc, argv); - - std::cout << "starting the exact Riemann solver..." << std::endl; - - // initialize the Castro runtime parameters + exact_riemann(); - Castro castro; - castro.doit(); } From a4a16d352c95772ed11fce64d28289df9a5007ec Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 12:48:45 -0400 Subject: [PATCH 04/25] this works now --- Util/exact_riemann/_prob_params | 2 -- Util/exact_riemann/exact_riemann.H | 4 +++- Util/exact_riemann/inputs.stellarcollapse | 4 ++-- Util/exact_riemann/inputs.test1.helm | 1 + Util/exact_riemann/riemann_rarefaction.H | 8 ++++---- Util/exact_riemann/riemann_sample.H | 6 +++--- Util/exact_riemann/riemann_shock.H | 4 ++-- Util/exact_riemann/riemann_star_state.H | 4 ++-- 8 files changed, 17 insertions(+), 16 deletions(-) diff --git a/Util/exact_riemann/_prob_params b/Util/exact_riemann/_prob_params index 9c43f0d3b5..b8cf88be92 100644 --- a/Util/exact_riemann/_prob_params +++ b/Util/exact_riemann/_prob_params @@ -19,5 +19,3 @@ t real 0.2d0 y npts integer 128 y use_Tinit integer 0 y - -initial_temp_guess real 1.0d5 y diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index 39deafe8d4..52b44db51a 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -5,6 +5,8 @@ #include +#include + #include #include #include @@ -88,7 +90,7 @@ exact_riemann() { for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn_s[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); diff --git a/Util/exact_riemann/inputs.stellarcollapse b/Util/exact_riemann/inputs.stellarcollapse index 92dd346435..412b02c272 100644 --- a/Util/exact_riemann/inputs.stellarcollapse +++ b/Util/exact_riemann/inputs.stellarcollapse @@ -8,8 +8,6 @@ problem.T_r = 5.e10 problem.use_Tinit = 1 -problem.initial_temp_guess = 1.e5 - problem.riemann_max_iter = 20 problem.xmin = 0.0 @@ -20,5 +18,7 @@ problem.t = 1.e-6 problem.npts = 128 +castro.T_guess = 1.e5 + eos.eos_file = LS220_240r_140t_50y_analmu_20120628_SVNr28.h5 diff --git a/Util/exact_riemann/inputs.test1.helm b/Util/exact_riemann/inputs.test1.helm index db1db7d1ce..b36d2f0bab 100644 --- a/Util/exact_riemann/inputs.test1.helm +++ b/Util/exact_riemann/inputs.test1.helm @@ -18,3 +18,4 @@ problem.t = 0.0008 problem.npts = 128 +castro.T_guess = 1.e5 diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H index c5d98c6f87..fd2e266342 100644 --- a/Util/exact_riemann/riemann_rarefaction.H +++ b/Util/exact_riemann/riemann_rarefaction.H @@ -185,7 +185,7 @@ rarefaction(const amrex::Real pstar, // initial guess at integration step amrex::Real dp = (pstar - p_s) / static_cast(100); - amrex::Real T = problem::initial_temp_guess; + amrex::Real T = castro::T_guess; // adaptive RK4 loop @@ -253,7 +253,7 @@ rarefaction(const amrex::Real pstar, for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -311,7 +311,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -389,7 +389,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index 9db07b34df..8f0350fa4c 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -25,7 +25,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn_l[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -36,7 +36,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn_r[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -129,7 +129,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index c91ac466de..3a74b6e3d3 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -26,7 +26,7 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar, for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -113,7 +113,7 @@ shock(const amrex::Real pstar, for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 67eeabada7..2448247129 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -28,7 +28,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn_l[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); @@ -42,7 +42,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: for (int n = 0; n < NumSpec; ++n) { eos_state.xn[n] = xn_r[n]; } - eos_state.T = problem::initial_temp_guess; + eos_state.T = castro::T_guess; eos(eos_input_rp, eos_state); From b0b48bbd6a0aeb69e6d482933307218b2df06656 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 12:51:41 -0400 Subject: [PATCH 05/25] add missing file --- Util/exact_riemann/struct_params.H | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 Util/exact_riemann/struct_params.H diff --git a/Util/exact_riemann/struct_params.H b/Util/exact_riemann/struct_params.H new file mode 100644 index 0000000000..8ad1776283 --- /dev/null +++ b/Util/exact_riemann/struct_params.H @@ -0,0 +1,8 @@ +#ifndef STRUCT_PARAMS_H +#define STRUCT_PARAMS_H + +#include + +inline params_t params; + +#endif From 8b810445d5b1a6158197f1fc14cee0257bb21126 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 15:27:02 -0400 Subject: [PATCH 06/25] unused header --- Util/exact_riemann/exact_riemann.H | 1 - 1 file changed, 1 deletion(-) diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index 52b44db51a..e99bf69a4d 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -9,7 +9,6 @@ #include #include -#include using namespace amrex::literals; From 372229c9eab929175e8b74c9a71895ade27e1349 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 16:44:23 -0400 Subject: [PATCH 07/25] rename some Riemann runtime parameters now they have "riemann" in their name, and also some drop "cg" so that it is clear that they can apply to other solvers --- Source/driver/Castro.cpp | 8 ++++---- Source/driver/_cpp_parameters | 15 +++++++------- Source/hydro/riemann_2shock_solvers.H | 28 +++++++++++++-------------- 3 files changed, 25 insertions(+), 26 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index d79f4d4935..f478f8402d 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -392,12 +392,12 @@ Castro::read_params () #endif if (riemann_solver == 1) { - if (cg_maxiter > HISTORY_SIZE) { - amrex::Error("error in riemanncg: cg_maxiter > HISTORY_SIZE"); + if (riemann_shock_maxiter > HISTORY_SIZE) { + amrex::Error("riemann_shock_maxiter > HISTORY_SIZE"); } - if (cg_blend == 2 && cg_maxiter < 5) { - amrex::Error("Error: need cg_maxiter >= 5 to do a bisection search on secant iteration failure."); + if (riemann_cg_blend == 2 && riemann_shock_maxiter < 5) { + amrex::Error("Error: need riemann_shock_maxiter >= 5 to do a bisection search on secant iteration failure."); } } #endif diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index a342478def..f5cdc29bd3 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -114,20 +114,19 @@ hybrid_riemann bool 0 # 2: HLLC riemann_solver int 0 -# for the Colella \& Glaz Riemann solver, the maximum number -# of iterations to take when solving for the star state -cg_maxiter int 12 +# maximum number of iterations to used in the Riemann solver +# when solving for the star state +riemann_shock_maxiter int 12 -# for the Colella \& Glaz Riemann solver, the tolerance to -# demand in finding the star state -cg_tol Real 1.0e-5 +# tolerance to use when finding the star stat +riemann_pstar_tol Real 1.0e-5 # for the Colella \& Glaz Riemann solver, what to do if # we do not converge to a solution for the star state. # 0 = do nothing; print iterations and exit # 1 = revert to the original guess for p-star -# 2 = do a bisection search for another 2 * cg_maxiter iterations. -cg_blend int 2 +# 2 = do a bisection search for another 2 * riemann_shock_maxiter iterations. +riemann_cg_blend int 2 # flatten the reconstructed profiles around shocks to prevent them # from becoming too thin diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H index 208e9878e0..f81b495f29 100644 --- a/Source/hydro/riemann_2shock_solvers.H +++ b/Source/hydro/riemann_2shock_solvers.H @@ -58,7 +58,7 @@ namespace TwoShock { const amrex::Real ur, const amrex::Real pr, const amrex::Real taur, const amrex::Real gamer, const amrex::Real clsqr, const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax, - const int lcg_maxiter, const amrex::Real lcg_tol, + const int lriemann_shock_maxiter, const amrex::Real lriemann_pstar_tol, amrex::Real& pstar, amrex::Real& gamstar, bool& converged, amrex::GpuArray& pstar_hist_extra) { @@ -106,7 +106,7 @@ namespace TwoShock { converged = false; amrex::Real pstar_c = 0.0; - for (int iter = 0; iter < PSTAR_BISECT_FACTOR*lcg_maxiter; iter++) { + for (int iter = 0; iter < PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) { pstar_c = 0.5_rt * (pstar_lo + pstar_hi); pstar_hist_extra[iter] = pstar_c; @@ -125,7 +125,7 @@ namespace TwoShock { amrex::Real f_c = ustar_l - ustar_r; - if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lcg_tol * pstar_c ) { + if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lriemann_pstar_tol * pstar_c ) { converged = true; break; } @@ -239,7 +239,7 @@ namespace TwoShock { bool converged = false; int iter = 0; - while ((iter < cg_maxiter && !converged) || iter < 2) { + while ((iter < riemann_shock_maxiter && !converged) || iter < 2) { wsqge(ql.p, taul, gamel, gdot, gamstar, gmin, gmax, clsql, pstar, wlsq); @@ -281,7 +281,7 @@ namespace TwoShock { pstar = amrex::max(pstar, small_pres); amrex::Real err = std::abs(pstar - pstar_old); - if (err < cg_tol*pstar) { + if (err < riemann_pstar_tol * pstar) { converged = true; } @@ -300,11 +300,11 @@ namespace TwoShock { if (!converged) { - if (cg_blend == 0) { + if (riemann_cg_blend == 0) { #ifndef AMREX_USE_GPU std::cout << "pstar history: " << std::endl; - for (int iter_l=0; iter_l < cg_maxiter; iter_l++) { + for (int iter_l=0; iter_l < riemann_shock_maxiter; iter_l++) { std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; } @@ -316,11 +316,11 @@ namespace TwoShock { amrex::Error("ERROR: non-convergence in the Riemann solver"); #endif - } else if (cg_blend == 1) { + } else if (riemann_cg_blend == 1) { pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); - } else if (cg_blend == 2) { + } else if (riemann_cg_blend == 2) { // we don't store the history if we are in CUDA, so // we can't do this @@ -328,7 +328,7 @@ namespace TwoShock { // first try to find a reasonable bounds amrex::Real pstarl = 1.e200; amrex::Real pstaru = -1.e200; - for (int n = cg_maxiter-6; n < cg_maxiter; n++) { + for (int n = riemann_shock_maxiter-6; n < riemann_shock_maxiter; n++) { pstarl = amrex::min(pstarl, pstar_hist[n]); pstaru = amrex::max(pstaru, pstar_hist[n]); } @@ -342,17 +342,17 @@ namespace TwoShock { ql.un, ql.p, taul, gamel, clsql, qr.un, qr.p, taur, gamer, clsqr, gdot, gmin, gmax, - cg_maxiter, cg_tol, + riemann_shock_maxiter, riemann_pstar_tol, pstar, gamstar, converged, pstar_hist_extra); if (!converged) { std::cout << "pstar history: " << std::endl; - for (int iter_l = 0; iter_l < cg_maxiter; iter_l++) { + for (int iter_l = 0; iter_l < riemann_shock_maxiter; iter_l++) { std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; } std::cout << "pstar extra history: " << std::endl; - for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR*cg_maxiter; iter_l++) { + for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) { std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; } @@ -368,7 +368,7 @@ namespace TwoShock { } else { #ifndef AMREX_USE_GPU - amrex::Error("ERROR: unrecognized cg_blend option."); + amrex::Error("ERROR: unrecognized riemann_cg_blend option."); #endif } From 1f3a733adc6c47f028b9ba48733d74b314e04df1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 16:51:47 -0400 Subject: [PATCH 08/25] fix runtime params --- Exec/science/flame_wave/ci-benchmarks/job_info_params.txt | 6 +++--- Exec/science/wdmerger/inputs | 2 +- Exec/science/wdmerger/tests/he_double_det/inputs_pakmor | 2 +- .../wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc | 2 +- .../wdmerger/tests/he_double_det/inputs_pakmor_strang | 2 +- Exec/science/wdmerger/tests/he_double_det/inputs_scaling | 2 +- .../wdmerger/tests/wdmerger_collision/inputs_2d_collision | 2 +- Exec/science/wdmerger/tests/wdmerger_collision_1D/inputs | 2 +- .../tests/wdmerger_retry/inputs_test_wdmerger_retry | 2 +- 9 files changed, 11 insertions(+), 11 deletions(-) diff --git a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt index 127b5bcb64..252b98a482 100644 --- a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -54,9 +54,9 @@ castro.plm_limiter = 2 castro.hybrid_riemann = 0 castro.riemann_solver = 0 - castro.cg_maxiter = 12 - castro.cg_tol = 1e-05 - castro.cg_blend = 2 + castro.riemann_shock_maxiter = 12 + castro.riemann_pstar_tol = 1e-05 + castro.riemann_cg_blend = 2 castro.use_flattening = 1 castro.transverse_use_eos = 0 castro.transverse_reset_density = 1 diff --git a/Exec/science/wdmerger/inputs b/Exec/science/wdmerger/inputs index 5def509a0d..f9b22104b2 100644 --- a/Exec/science/wdmerger/inputs +++ b/Exec/science/wdmerger/inputs @@ -188,7 +188,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor index 2ff2712602..a7028d1544 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor @@ -200,7 +200,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc index db97ad3dac..6316753695 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc @@ -206,7 +206,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang index a57a028d0d..b56d9ba2a5 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang @@ -206,7 +206,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_scaling b/Exec/science/wdmerger/tests/he_double_det/inputs_scaling index 03c95c1636..990c81fd2a 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_scaling +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_scaling @@ -203,7 +203,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/wdmerger_collision/inputs_2d_collision b/Exec/science/wdmerger/tests/wdmerger_collision/inputs_2d_collision index 8a2bca0293..60a7db22d6 100644 --- a/Exec/science/wdmerger/tests/wdmerger_collision/inputs_2d_collision +++ b/Exec/science/wdmerger/tests/wdmerger_collision/inputs_2d_collision @@ -153,7 +153,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/wdmerger_collision_1D/inputs b/Exec/science/wdmerger/tests/wdmerger_collision_1D/inputs index 6dcabf249c..9f40085682 100644 --- a/Exec/science/wdmerger/tests/wdmerger_collision_1D/inputs +++ b/Exec/science/wdmerger/tests/wdmerger_collision_1D/inputs @@ -142,7 +142,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 diff --git a/Exec/science/wdmerger/tests/wdmerger_retry/inputs_test_wdmerger_retry b/Exec/science/wdmerger/tests/wdmerger_retry/inputs_test_wdmerger_retry index 0d63ab606c..a6eef171b3 100644 --- a/Exec/science/wdmerger/tests/wdmerger_retry/inputs_test_wdmerger_retry +++ b/Exec/science/wdmerger/tests/wdmerger_retry/inputs_test_wdmerger_retry @@ -133,7 +133,7 @@ castro.riemann_solver = 0 # For the CG Riemann solver, we need to tell the solver not to quit when # the iterations don't converge, but instead to do additional bisection iteration. -castro.cg_blend = 2 +castro.riemann_cg_blend = 2 # Use a lagged predictor estimate of the source terms in the hydro castro.source_term_predictor = 1 From f6363f2240c17af13449dd57feab38fe79f2d9e6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 04:55:24 -0400 Subject: [PATCH 09/25] work on using Source/hydro riemann parameters --- Util/exact_riemann/main.cpp | 4 +--- Util/exact_riemann/riemann_shock.H | 9 +++++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp index 7f75f6fc36..4707027763 100644 --- a/Util/exact_riemann/main.cpp +++ b/Util/exact_riemann/main.cpp @@ -41,9 +41,7 @@ int main(int argc, char *argv[]) { network_init(); #endif - amrex::Real small_temp = 1.e-200; - amrex::Real small_dens = 1.e-200; - eos_init(small_temp, small_dens); + eos_init(castro::small_temp, castro::small_dens); exact_riemann(); diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index 3a74b6e3d3..8a20c598a0 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -6,6 +6,7 @@ #include #include +#include AMREX_INLINE void @@ -62,7 +63,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar, converged = false; int iter = 1; - while (! converged && iter < riemann_solver::max_iters) { + while (! converged && iter < castro::riemann_shock_maxiter) { amrex::Real rhostar_s, f, fprime; @@ -71,7 +72,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar, amrex::Real dW = -f / fprime; - if (std::abs(dW) < riemann_solver::tol * W_s) { + if (std::abs(dW) < castro::riemann_pstar_tol * W_s) { converged = true; } @@ -147,7 +148,7 @@ shock(const amrex::Real pstar, amrex::Real rhostar_s; if (taustar_s < 0.0_rt) { - rhostar_s = riemann_solver::smallrho; + rhostar_s = riemann_constants::small; W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s)); } @@ -162,7 +163,7 @@ shock(const amrex::Real pstar, // now did we converge? if (! converged) { - for (int i = 0; i < riemann_solver::max_iters; ++i) { + for (int i = 0; i < castro::riemann_riemann_solver::max_iters; ++i) { std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl; } From 03b9653a2198a9567501e2d2d83a1d15b3649020 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 05:15:17 -0400 Subject: [PATCH 10/25] move some Riemann history constexpr to riemann_constants namespace --- Source/driver/Castro.H | 6 ++---- Source/driver/Castro.cpp | 2 +- Source/hydro/riemann_2shock_solvers.H | 10 +++++----- Source/hydro/riemann_type.H | 10 +++++++--- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 30976c248d..f5447c0427 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -16,10 +16,6 @@ using namespace castro; #include -constexpr int HISTORY_SIZE=40; -constexpr int PSTAR_BISECT_FACTOR = 5; - - #include #include #include @@ -50,6 +46,8 @@ constexpr int PSTAR_BISECT_FACTOR = 5; #include +#include + using std::istream; using std::ostream; diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index f478f8402d..adc4d97fe9 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -392,7 +392,7 @@ Castro::read_params () #endif if (riemann_solver == 1) { - if (riemann_shock_maxiter > HISTORY_SIZE) { + if (riemann_shock_maxiter > riemann_constants::HISTORY_SIZE) { amrex::Error("riemann_shock_maxiter > HISTORY_SIZE"); } diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H index f81b495f29..9fc55d5839 100644 --- a/Source/hydro/riemann_2shock_solvers.H +++ b/Source/hydro/riemann_2shock_solvers.H @@ -60,7 +60,7 @@ namespace TwoShock { const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax, const int lriemann_shock_maxiter, const amrex::Real lriemann_pstar_tol, amrex::Real& pstar, amrex::Real& gamstar, - bool& converged, amrex::GpuArray& pstar_hist_extra) { + bool& converged, amrex::GpuArray& pstar_hist_extra) { // we want to zero // f(p*) = u*_l(p*) - u*_r(p*) @@ -106,7 +106,7 @@ namespace TwoShock { converged = false; amrex::Real pstar_c = 0.0; - for (int iter = 0; iter < PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) { + for (int iter = 0; iter < riemann_constants::PSTAR_BISECT_FACTOR * lriemann_shock_maxiter; iter++) { pstar_c = 0.5_rt * (pstar_lo + pstar_hi); pstar_hist_extra[iter] = pstar_c; @@ -167,7 +167,7 @@ namespace TwoShock { constexpr amrex::Real weakwv = 1.e-3_rt; #ifndef AMREX_USE_GPU - amrex::GpuArray pstar_hist; + amrex::GpuArray pstar_hist; #endif @@ -336,7 +336,7 @@ namespace TwoShock { pstarl = amrex::max(pstarl, small_pres); pstaru = amrex::max(pstaru, small_pres); - amrex::GpuArray pstar_hist_extra; + amrex::GpuArray pstar_hist_extra; pstar_bisection(pstarl, pstaru, ql.un, ql.p, taul, gamel, clsql, @@ -352,7 +352,7 @@ namespace TwoShock { std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; } std::cout << "pstar extra history: " << std::endl; - for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) { + for (int iter_l = 0; iter_l < riemann_constants::PSTAR_BISECT_FACTOR * riemann_shock_maxiter; iter_l++) { std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; } diff --git a/Source/hydro/riemann_type.H b/Source/hydro/riemann_type.H index c8be48f47a..5ca1a9778b 100644 --- a/Source/hydro/riemann_type.H +++ b/Source/hydro/riemann_type.H @@ -6,9 +6,13 @@ using namespace amrex::literals; namespace riemann_constants { - const amrex::Real smlp1 = 1.e-10_rt; - const amrex::Real small = 1.e-8_rt; - const amrex::Real smallu = 1.e-12_rt; + constexpr amrex::Real smlp1 = 1.e-10_rt; + constexpr amrex::Real small = 1.e-8_rt; + constexpr amrex::Real smallu = 1.e-12_rt; + constexpr int HISTORY_SIZE=40; + constexpr int PSTAR_BISECT_FACTOR = 5; + + } From 9e7e820e74a8a2376a4e9ff29e9fee2297059900 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 05:40:52 -0400 Subject: [PATCH 11/25] split constants into riemann_constants.H allow the exact Riemann solver to use the new constants --- Source/driver/Castro.cpp | 2 ++ Source/hydro/Make.package | 1 + Source/hydro/riemann_2shock_solvers.H | 1 + Source/hydro/riemann_constants.H | 16 ++++++++++++++++ Source/hydro/riemann_type.H | 11 ----------- Util/exact_riemann/GNUmakefile | 2 +- Util/exact_riemann/riemann_shock.H | 6 +++--- 7 files changed, 24 insertions(+), 15 deletions(-) create mode 100644 Source/hydro/riemann_constants.H diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index adc4d97fe9..e79cf0574e 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -58,6 +58,8 @@ #include #include +#include + using namespace amrex; bool Castro::signalStopJob = false; diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package index 6b6ba6c7e8..c5813e27c3 100644 --- a/Source/hydro/Make.package +++ b/Source/hydro/Make.package @@ -30,6 +30,7 @@ CEXE_sources += riemann.cpp CEXE_headers += HLL_solvers.H CEXE_headers += riemann_solvers.H CEXE_headers += riemann_2shock_solvers.H +CEXE_headers += riemann_constants.H CEXE_headers += riemann_type.H CEXE_headers += slope.H CEXE_headers += reconstruction.H diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H index 9fc55d5839..6a70d0e117 100644 --- a/Source/hydro/riemann_2shock_solvers.H +++ b/Source/hydro/riemann_2shock_solvers.H @@ -5,6 +5,7 @@ using namespace amrex::literals; +#include #include #ifdef RADIATION #include diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H new file mode 100644 index 0000000000..65ca6641c8 --- /dev/null +++ b/Source/hydro/riemann_constants.H @@ -0,0 +1,16 @@ +#ifndef RIEMANN_CONSTANTS +#define RIEMANN_CONSTANTS + +#include + +using namespace amrex::literals; + +namespace riemann_constants { + constexpr amrex::Real smlp1 = 1.e-10_rt; + constexpr amrex::Real small = 1.e-8_rt; + constexpr amrex::Real smallu = 1.e-12_rt; + constexpr int HISTORY_SIZE=40; + constexpr int PSTAR_BISECT_FACTOR = 5; +} + +#endif diff --git a/Source/hydro/riemann_type.H b/Source/hydro/riemann_type.H index 5ca1a9778b..658befe151 100644 --- a/Source/hydro/riemann_type.H +++ b/Source/hydro/riemann_type.H @@ -5,17 +5,6 @@ using namespace amrex::literals; -namespace riemann_constants { - constexpr amrex::Real smlp1 = 1.e-10_rt; - constexpr amrex::Real small = 1.e-8_rt; - constexpr amrex::Real smallu = 1.e-12_rt; - constexpr int HISTORY_SIZE=40; - constexpr int PSTAR_BISECT_FACTOR = 5; - - -} - - struct RiemannState { amrex::Real rho; diff --git a/Util/exact_riemann/GNUmakefile b/Util/exact_riemann/GNUmakefile index 8be2709e1a..fbf934926a 100644 --- a/Util/exact_riemann/GNUmakefile +++ b/Util/exact_riemann/GNUmakefile @@ -33,6 +33,6 @@ Blocs := . # we explicitly want runtime_parameters.H so we have access to # Castro's runtime parameter system -Blocs += $(CASTRO_HOME)/Source/driver +Blocs += $(CASTRO_HOME)/Source/driver $(CASTRO_HOME)/Source/hydro include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index 8a20c598a0..a6bb04eb9e 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -6,7 +6,7 @@ #include #include -#include +#include AMREX_INLINE void @@ -97,7 +97,7 @@ shock(const amrex::Real pstar, amrex::ignore_unused(u_s); - amrex::Real rhostar_hist[riemann_solver::max_iters], Ws_hist[riemann_solver::max_iters]; + amrex::Real rhostar_hist[riemann_constants::HISTORY_SIZE], Ws_hist[riemann_constants::HISTORY_SIZE]; // compute the Z_s function for a shock following C&G Eq. 20 and // 23. Here the "_s" variables are the state (L or R) that we are @@ -163,7 +163,7 @@ shock(const amrex::Real pstar, // now did we converge? if (! converged) { - for (int i = 0; i < castro::riemann_riemann_solver::max_iters; ++i) { + for (int i = 0; i < castro::riemann_shock_maxiter; ++i) { std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl; } From 9b1cf4d47e6457ca459b124dc6d07d25843bbb98 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 07:05:41 -0400 Subject: [PATCH 12/25] remove unneeded header --- Source/driver/Castro.H | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index f5447c0427..8c7b9c13f8 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -46,8 +46,6 @@ using namespace castro; #include -#include - using std::istream; using std::ostream; From 0a6ce88d0ec1f8fe4dabba82c5c0c35e11900083 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 07:56:06 -0400 Subject: [PATCH 13/25] more work on getting the Riemann solver to fully use Castro stuff --- Source/hydro/riemann_constants.H | 3 +++ Util/exact_riemann/Make.package | 1 - Util/exact_riemann/inputs.test1.helm | 1 + Util/exact_riemann/inputs.test2.helm | 2 ++ Util/exact_riemann/riemann_rarefaction.H | 14 +++++--------- Util/exact_riemann/riemann_sample.H | 2 +- Util/exact_riemann/riemann_shock.H | 2 +- Util/exact_riemann/riemann_star_state.H | 11 ++++++----- Util/exact_riemann/riemann_support.H | 14 -------------- 9 files changed, 19 insertions(+), 31 deletions(-) delete mode 100644 Util/exact_riemann/riemann_support.H diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H index 65ca6641c8..ad8fa2a141 100644 --- a/Source/hydro/riemann_constants.H +++ b/Source/hydro/riemann_constants.H @@ -9,6 +9,9 @@ namespace riemann_constants { constexpr amrex::Real smlp1 = 1.e-10_rt; constexpr amrex::Real small = 1.e-8_rt; constexpr amrex::Real smallu = 1.e-12_rt; + constexpr amrex::Real riemann_integral_tol = 1.e-8_rt; + constexpr amrex::Real riemann_u_tol = 1.e-6_rt; + constexpr amrex::Real riemann_p_tol = 1.e-8_rt; constexpr int HISTORY_SIZE=40; constexpr int PSTAR_BISECT_FACTOR = 5; } diff --git a/Util/exact_riemann/Make.package b/Util/exact_riemann/Make.package index 2cc1befab7..4916ffc34c 100644 --- a/Util/exact_riemann/Make.package +++ b/Util/exact_riemann/Make.package @@ -3,7 +3,6 @@ CEXE_sources += main.cpp CEXE_headers += exact_riemann.H CEXE_headers += riemann_star_state.H CEXE_headers += riemann_sample.H -CEXE_headers += riemann_support.H CEXE_sources += extern_parameters.cpp CEXE_headers += extern_parameters.H diff --git a/Util/exact_riemann/inputs.test1.helm b/Util/exact_riemann/inputs.test1.helm index b36d2f0bab..2810e90046 100644 --- a/Util/exact_riemann/inputs.test1.helm +++ b/Util/exact_riemann/inputs.test1.helm @@ -19,3 +19,4 @@ problem.t = 0.0008 problem.npts = 128 castro.T_guess = 1.e5 +castro.riemann_shock_maxiter = 20 \ No newline at end of file diff --git a/Util/exact_riemann/inputs.test2.helm b/Util/exact_riemann/inputs.test2.helm index 4a96e631a3..79bcd27016 100644 --- a/Util/exact_riemann/inputs.test2.helm +++ b/Util/exact_riemann/inputs.test2.helm @@ -15,3 +15,5 @@ problem.xjump = 0.5e5 problem.t = 0.00008 problem.npts = 128 + +castro.riemann_shock_maxiter = 100 diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H index fd2e266342..8cdefed510 100644 --- a/Util/exact_riemann/riemann_rarefaction.H +++ b/Util/exact_riemann/riemann_rarefaction.H @@ -174,8 +174,6 @@ rarefaction(const amrex::Real pstar, const double S1{0.9}; const double S2{4.0}; - constexpr amrex::Real tol = 1.e-5; - // initial conditions amrex::Real tau = 1.0_rt / rho_s; @@ -206,7 +204,7 @@ rarefaction(const amrex::Real pstar, // take a step - while (rel_error > tol) { + while (rel_error > riemann_constants::riemann_integral_tol) { dp = dp_new; if (p + dp < pstar) { dp = pstar - p; @@ -232,7 +230,7 @@ rarefaction(const amrex::Real pstar, // adaptive step estimate - double dp_est = S1 * dp * std::pow(std::abs(tol/rel_error), 0.2); + double dp_est = S1 * dp * std::pow(std::abs(riemann_constants::riemann_integral_tol/rel_error), 0.2); dp_new = dp_est; //std::clamp(S1*dp_est, dp/S2, S2*dp); } @@ -295,8 +293,6 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re const double S1{0.9}; - constexpr amrex::Real tol = 1.e-8; - // initial conditions amrex::Real tau = 1.0_rt / rho_s; @@ -343,7 +339,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re // take a step - while (rel_error > tol) { + while (rel_error > riemann_constants::riemann_integral_tol) { du = du_new; // take 2 half steps @@ -368,7 +364,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re amrex::Real du_est{}; if (rel_error > 0) { - du_est = S1 * du * std::pow(std::abs(tol / rel_error), 0.2); + du_est = S1 * du * std::pow(std::abs(riemann_constants::riemann_integral_tol / rel_error), 0.2); } else { du_est = 10.0 * du; } @@ -417,7 +413,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re } } - if (std::abs(du_new) < riemann_solver::tol * std::abs(u) + riemann_solver::tol) { + if (std::abs(du_new) < riemann_constants::riemann_u_tol * std::abs(u) + riemann_constants::riemann_u_tol) { finished = true; } diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index 8f0350fa4c..d123fd046b 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -3,7 +3,7 @@ #include -#include +#include #include using namespace amrex::literals; diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index a6bb04eb9e..ff21ea5932 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -134,7 +134,7 @@ shock(const amrex::Real pstar, // just doesn't work. In this case, we use the ideas from CG, Eq. 35, and // take W = sqrt(gamma p rho) - if (pstar - p_s < riemann_solver::tol_p * p_s) { + if (pstar - p_s < riemann_constants::riemann_p_tol * p_s) { W_s = std::sqrt(eos_state.gam1 * p_s * rho_s); } else { W_s = std::sqrt((pstar - p_s) * diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 2448247129..2eb83221d5 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -3,7 +3,7 @@ #include -#include +#include #include #include @@ -108,7 +108,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: bool converged = false; int iter = 1; - while (! converged && iter < riemann_solver::max_iters) { + while (! converged && iter < castro::riemann_shock_maxiter) { // compute Z_l, W_l and Z_r, W_r -- the form of these depend // on whether the wave is a shock or a rarefaction @@ -146,13 +146,14 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: amrex::Real err1 = std::abs(ustar_r - ustar_l); amrex::Real err2 = pstar_new - pstar; - if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && - err2 < tol * pstar) { + std::cout << "iter = " << iter << " " << err1 / std::max(std::abs(ustar_l), std::abs(ustar_r)) << " " << err2 / pstar << std::endl; + + if (err2 < tol * pstar) { converged = true; } // get ready for the next iteration - pstar = pstar_new; + pstar = std::clamp(pstar_new, 0.5 * pstar, 2.0 * pstar); iter++; } diff --git a/Util/exact_riemann/riemann_support.H b/Util/exact_riemann/riemann_support.H deleted file mode 100644 index 56fd6b8aaf..0000000000 --- a/Util/exact_riemann/riemann_support.H +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef RIEMANN_SUPPORT_H -#define RIEMANN_SUPPORT_H - -#include - -namespace riemann_solver { - const int max_iters = 100; - - const amrex::Real smallrho = 1.e-5_rt; - const amrex::Real tol = 1.e-6_rt; - const amrex::Real tol_p = 1.e-6_rt; -} - -#endif From 17431baebae284a60d0eaecde8d08e1446c25c32 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 08:28:19 -0400 Subject: [PATCH 14/25] for pstar convergence in exact_riemann use std::abs() --- Util/exact_riemann/riemann_star_state.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 2448247129..5aeee19bcf 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -144,7 +144,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: // estimate the error in the current star solution amrex::Real err1 = std::abs(ustar_r - ustar_l); - amrex::Real err2 = pstar_new - pstar; + amrex::Real err2 = std::abs(pstar_new - pstar); if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && err2 < tol * pstar) { From 65a818c83dbc7061f1f671bc659cd28e3606f670 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 10:22:05 -0400 Subject: [PATCH 15/25] remove the ustar test --- Util/exact_riemann/README.md | 11 +++++++++++ Util/exact_riemann/riemann_star_state.H | 9 ++++----- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/Util/exact_riemann/README.md b/Util/exact_riemann/README.md index ada8e281ec..013a81693a 100644 --- a/Util/exact_riemann/README.md +++ b/Util/exact_riemann/README.md @@ -12,3 +12,14 @@ https://ui.adsabs.harvard.edu/abs/2015ApJS..216...31Z/abstract and more details are given there. To build the solver, simply type 'make' in this directory. + + +Notes: + +* ``test2`` is a vacuum + + The velocity in the star region should be 0, but it seems that + slight differences in the rarefaction integration from the left and + from the right give different roundoff, resulting in a small, finite + velocity. This can be controlled a bit by tightening the tolerance + in the rarefaction integration. diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 5aeee19bcf..31e1952722 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -142,12 +142,11 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: amrex::Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r); - // estimate the error in the current star solution - amrex::Real err1 = std::abs(ustar_r - ustar_l); - amrex::Real err2 = std::abs(pstar_new - pstar); + // estimate the error in the current pstar solution - if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && - err2 < tol * pstar) { + amrex::Real error = std::abs(pstar_new - pstar); + + if (error < tol * pstar) { converged = true; } From 22a9572ae78ab556fe0cd898122d35bd5e8babb1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 8 Jul 2024 08:49:02 -0400 Subject: [PATCH 16/25] fix header guard --- Source/hydro/riemann_constants.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H index 65ca6641c8..f06776238f 100644 --- a/Source/hydro/riemann_constants.H +++ b/Source/hydro/riemann_constants.H @@ -1,5 +1,5 @@ -#ifndef RIEMANN_CONSTANTS -#define RIEMANN_CONSTANTS +#ifndef RIEMANN_CONSTANTS_H +#define RIEMANN_CONSTANTS_H #include From 3fb9b516b8672a5fd0612ef6cdc7ab608bcf1caa Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 9 Jul 2024 09:58:40 -0400 Subject: [PATCH 17/25] more cleaning --- .../ci-benchmarks/test1-riemann.out | 90 +++++++++---------- Util/exact_riemann/exact_riemann.H | 9 +- Util/exact_riemann/inputs.bug8 | 1 + Util/exact_riemann/inputs.test1.helm | 1 - Util/exact_riemann/inputs.test2.helm | 2 - Util/exact_riemann/riemann_shock.H | 17 ++-- Util/exact_riemann/riemann_star_state.H | 7 +- 7 files changed, 64 insertions(+), 63 deletions(-) diff --git a/Util/exact_riemann/ci-benchmarks/test1-riemann.out b/Util/exact_riemann/ci-benchmarks/test1-riemann.out index 1624fcf02c..1c68b0fe98 100644 --- a/Util/exact_riemann/ci-benchmarks/test1-riemann.out +++ b/Util/exact_riemann/ci-benchmarks/test1-riemann.out @@ -68,51 +68,51 @@ 67 519531.25 3782769.06389 306131662.364 1.99799647911e+23 58852893.6708 9.19281449819e+16 1.50260347548 68 527343.75 3677627.42215 314046760.692 1.9150955856e+23 57947741.148 9.04497857546e+16 1.5041134824 69 535156.25 3573422.36381 322064216.687 1.83402150589e+23 57038358.4286 8.89637619826e+16 1.50565406349 - 70 542968.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 71 550781.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 72 558593.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 73 566406.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 74 574218.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 75 582031.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 76 589843.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 77 597656.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 78 605468.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 79 613281.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 80 621093.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 81 628906.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 82 636718.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 83 644531.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 84 652343.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 85 660156.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 86 667968.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 87 675781.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 88 683593.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 89 691406.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 90 699218.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 91 707031.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 92 714843.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 93 722656.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 94 730468.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 95 738281.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 96 746093.75 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 97 753906.25 3566771.39183 322581533.75 1.82889264698e+23 57003293.328 8.88686252995e+16 1.50575410519 - 98 761718.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 99 769531.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 100 777343.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 101 785156.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 102 792968.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 103 800781.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 104 808593.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 105 816406.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 106 824218.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 107 832031.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 108 839843.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 109 847656.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 110 855468.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 111 863281.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 112 871093.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 113 878906.25 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 - 114 886718.75 2952884.63392 322581533.75 1.82889264698e+23 758658822.183 1.09750124564e+17 1.49554976542 + 70 542968.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 71 550781.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 72 558593.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 73 566406.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 74 574218.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 75 582031.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 76 589843.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 77 597656.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 78 605468.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 79 613281.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 80 621093.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 81 628906.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 82 636718.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 83 644531.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 84 652343.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 85 660156.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 86 667968.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 87 675781.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 88 683593.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 89 691406.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 90 699218.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 91 707031.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 92 714843.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 93 722656.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 94 730468.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 95 738281.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 96 746093.75 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 97 753906.25 3566780.0111 322581309.173 1.82889079023e+23 56979962.0431 8.88683155688e+16 1.50575384117 + 98 761718.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 99 769531.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 100 777343.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 101 785156.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 102 792968.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 103 800781.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 104 808593.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 105 816406.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 106 824218.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 107 832031.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 108 839843.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 109 847656.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 110 855468.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 111 863281.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 112 871093.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 113 878906.25 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 + 114 886718.75 2952883.40957 322581309.173 1.82889079023e+23 758657908.288 1.09750048534e+17 1.49554982528 115 894531.25 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 116 902343.75 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 117 910156.25 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index e99bf69a4d..82dec000ad 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -50,10 +50,13 @@ exact_riemann() { amrex::Real ustar, pstar, W_l, W_r; - riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn, - problem::rho_r, problem::u_r, problem::p_r, xn, - ustar, pstar, W_l, W_r); + std::cout << "solving for star state: " << std::endl; + + auto niter = riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn, + problem::rho_r, problem::u_r, problem::p_r, xn, + ustar, pstar, W_l, W_r); + std::cout << "found star state, niter = " << niter << "; pstar, ustar = " << pstar << " " << ustar << std::endl; // find the solution as a function of xi = x/t std::ofstream ofile; diff --git a/Util/exact_riemann/inputs.bug8 b/Util/exact_riemann/inputs.bug8 index 8cf2b3486b..793aafe012 100644 --- a/Util/exact_riemann/inputs.bug8 +++ b/Util/exact_riemann/inputs.bug8 @@ -16,3 +16,4 @@ problem.t = 0.0008 problem.npts = 512 +castro.riemann_shock_maxiter = 100 \ No newline at end of file diff --git a/Util/exact_riemann/inputs.test1.helm b/Util/exact_riemann/inputs.test1.helm index 2810e90046..b36d2f0bab 100644 --- a/Util/exact_riemann/inputs.test1.helm +++ b/Util/exact_riemann/inputs.test1.helm @@ -19,4 +19,3 @@ problem.t = 0.0008 problem.npts = 128 castro.T_guess = 1.e5 -castro.riemann_shock_maxiter = 20 \ No newline at end of file diff --git a/Util/exact_riemann/inputs.test2.helm b/Util/exact_riemann/inputs.test2.helm index 79bcd27016..4a96e631a3 100644 --- a/Util/exact_riemann/inputs.test2.helm +++ b/Util/exact_riemann/inputs.test2.helm @@ -15,5 +15,3 @@ problem.xjump = 0.5e5 problem.t = 0.00008 problem.npts = 128 - -castro.riemann_shock_maxiter = 100 diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index ff21ea5932..8514ba5808 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -46,12 +46,12 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar, AMREX_INLINE -void +bool newton_shock(amrex::Real& W_s, const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, const amrex::Real* xn, amrex::Real* rhostar_hist, amrex::Real* Ws_hist, - eos_rep_t& eos_state, bool& converged) { + eos_rep_t& eos_state) { // Newton iterations -- we are zeroing the energy R-H jump condition // W^2 [e] = 1/2 [p^2] @@ -60,7 +60,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar, // // and then compute f' - converged = false; + bool converged = false; int iter = 1; while (! converged && iter < castro::riemann_shock_maxiter) { @@ -76,7 +76,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar, converged = true; } - W_s = amrex::min(2.0_rt * W_s, amrex::max(0.5_rt * W_s, W_s + dW)); + W_s = std::clamp(W_s + dW, 0.5_rt * W_s, 2.0_rt * W_s); // store some history @@ -85,6 +85,8 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar, iter++; } + + return converged; } AMREX_INLINE @@ -154,10 +156,9 @@ shock(const amrex::Real pstar, // newton - bool converged; - newton_shock(W_s, pstar, rho_s, p_s, e_s, xn, - rhostar_hist, Ws_hist, - eos_state, converged); + auto converged = newton_shock(W_s, pstar, rho_s, p_s, e_s, xn, + rhostar_hist, Ws_hist, + eos_state); // now did we converge? diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index b12d148868..2bf7c12a15 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -10,7 +10,7 @@ using namespace amrex::literals; AMREX_INLINE -void +int riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real p_l, const amrex::Real* xn_l, const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r, amrex::Real& ustar, amrex::Real& pstar, amrex::Real& W_l, amrex::Real& W_r) { @@ -79,8 +79,6 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: // find the exact pstar and ustar - std::cout << "solving for star state: " << std::endl; - // this procedure follows directly from Colella & Glaz 1985, section 1 // the basic idea is that we want to find the pstar that satisfies: // @@ -162,6 +160,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: ustar = 0.5_rt * (ustar_l + ustar_r); - std::cout << "found pstar, ustar: " << pstar << " " << ustar << std::endl; + return iter; + } #endif From 51313d2dc40f52a34c7cd377db86ea109e2fa066 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 9 Jul 2024 10:01:58 -0400 Subject: [PATCH 18/25] finish --- Util/exact_riemann/riemann_shock.H | 4 ++-- Util/exact_riemann/riemann_star_state.H | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index 8514ba5808..5cf4e5efe0 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -40,7 +40,7 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar, amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; fprime = 2.0_rt * W_s * (eos_state.e - e_s) - - 2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s; + 2.0_rt * dedrho_p * (pstar - p_s) * amrex::Math::powi<2>(rhostar_s) / W_s; } @@ -186,7 +186,7 @@ shock(const amrex::Real pstar, amrex::Real p_e = eos_state.dpdT / eos_state.dedT; amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; - amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho; + amrex::Real p_tau = -amrex::Math::powi<2>(rhostar_s) * p_rho; amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / ((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s)); diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 2bf7c12a15..b432a61d82 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -74,7 +74,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: pstar = ((W_r * p_l + W_l * p_r) + W_l * W_r * (u_l - u_r)) / (W_l + W_r); } - pstar = amrex::max(pstar, smallp); + pstar = std::max(pstar, smallp); // find the exact pstar and ustar From 2a1a272edaabbcf3050dae4fd7b77774b76ca7ee Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 10:12:51 -0400 Subject: [PATCH 19/25] rename --- Util/exact_riemann/Make.package | 6 ++++-- Util/exact_riemann/exact_riemann.H | 4 ++-- .../{riemann_rarefaction.H => exact_riemann_rarefaction.H} | 0 .../{riemann_sample.H => exact_riemann_sample.H} | 2 +- .../{riemann_shock.H => exact_riemann_shock.H} | 0 .../{riemann_star_state.H => exact_riemann_star_state.H} | 4 ++-- 6 files changed, 9 insertions(+), 7 deletions(-) rename Util/exact_riemann/{riemann_rarefaction.H => exact_riemann_rarefaction.H} (100%) rename Util/exact_riemann/{riemann_sample.H => exact_riemann_sample.H} (99%) rename Util/exact_riemann/{riemann_shock.H => exact_riemann_shock.H} (100%) rename Util/exact_riemann/{riemann_star_state.H => exact_riemann_star_state.H} (98%) diff --git a/Util/exact_riemann/Make.package b/Util/exact_riemann/Make.package index 4916ffc34c..8e4f2a074e 100644 --- a/Util/exact_riemann/Make.package +++ b/Util/exact_riemann/Make.package @@ -1,8 +1,10 @@ CEXE_sources += main.cpp CEXE_headers += exact_riemann.H -CEXE_headers += riemann_star_state.H -CEXE_headers += riemann_sample.H +CEXE_headers += exact_riemann_star_state.H +CEXE_headers += exact_riemann_sample.H +CEXE_headers += exact_riemann_shock.H +CEXE_headers += exact_riemann_rarefaction.H CEXE_sources += extern_parameters.cpp CEXE_headers += extern_parameters.H diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index 82dec000ad..b1cf14e832 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -7,8 +7,8 @@ #include -#include -#include +#include +#include using namespace amrex::literals; diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/exact_riemann_rarefaction.H similarity index 100% rename from Util/exact_riemann/riemann_rarefaction.H rename to Util/exact_riemann/exact_riemann_rarefaction.H diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/exact_riemann_sample.H similarity index 99% rename from Util/exact_riemann/riemann_sample.H rename to Util/exact_riemann/exact_riemann_sample.H index d123fd046b..c255ea1cce 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/exact_riemann_sample.H @@ -4,7 +4,7 @@ #include #include -#include +#include using namespace amrex::literals; diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/exact_riemann_shock.H similarity index 100% rename from Util/exact_riemann/riemann_shock.H rename to Util/exact_riemann/exact_riemann_shock.H diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/exact_riemann_star_state.H similarity index 98% rename from Util/exact_riemann/riemann_star_state.H rename to Util/exact_riemann/exact_riemann_star_state.H index b432a61d82..7a421c1512 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/exact_riemann_star_state.H @@ -4,8 +4,8 @@ #include #include -#include -#include +#include +#include using namespace amrex::literals; From 6bfab826acd103f1cffdeea343839c9bc753f467 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 10:19:20 -0400 Subject: [PATCH 20/25] some moving --- .../hydro}/exact_riemann_rarefaction.H | 0 .../hydro}/exact_riemann_sample.H | 0 .../hydro}/exact_riemann_shock.H | 0 .../hydro}/exact_riemann_star_state.H | 0 Util/exact_riemann/README.md | 13 +++++++------ 5 files changed, 7 insertions(+), 6 deletions(-) rename {Util/exact_riemann => Source/hydro}/exact_riemann_rarefaction.H (100%) rename {Util/exact_riemann => Source/hydro}/exact_riemann_sample.H (100%) rename {Util/exact_riemann => Source/hydro}/exact_riemann_shock.H (100%) rename {Util/exact_riemann => Source/hydro}/exact_riemann_star_state.H (100%) diff --git a/Util/exact_riemann/exact_riemann_rarefaction.H b/Source/hydro/exact_riemann_rarefaction.H similarity index 100% rename from Util/exact_riemann/exact_riemann_rarefaction.H rename to Source/hydro/exact_riemann_rarefaction.H diff --git a/Util/exact_riemann/exact_riemann_sample.H b/Source/hydro/exact_riemann_sample.H similarity index 100% rename from Util/exact_riemann/exact_riemann_sample.H rename to Source/hydro/exact_riemann_sample.H diff --git a/Util/exact_riemann/exact_riemann_shock.H b/Source/hydro/exact_riemann_shock.H similarity index 100% rename from Util/exact_riemann/exact_riemann_shock.H rename to Source/hydro/exact_riemann_shock.H diff --git a/Util/exact_riemann/exact_riemann_star_state.H b/Source/hydro/exact_riemann_star_state.H similarity index 100% rename from Util/exact_riemann/exact_riemann_star_state.H rename to Source/hydro/exact_riemann_star_state.H diff --git a/Util/exact_riemann/README.md b/Util/exact_riemann/README.md index 013a81693a..77dbd4d069 100644 --- a/Util/exact_riemann/README.md +++ b/Util/exact_riemann/README.md @@ -1,9 +1,10 @@ -This is an exact Riemann solver for a general equation of state. It -follows the outline for Colella & Glaz 1985, section 1. This is too -slow to use in an actual hydro run, but instead is intended to -generate exact solutions to the Riemann problem for comparison with -Castro shocktube output. Several inputs files for Helmholtz EOS-based -shocktubes are provided. +This is a driver for the exact Riemann solver for a general equation +of state. The main implementation is in Source/hydro. + +The exact Riemann solver follows the outline for Colella & Glaz 1985, +section 1 and this driver is intended to generate exact solutions to +the Riemann problem for comparison with Castro shocktube output. +Several inputs files for Helmholtz EOS-based shocktubes are provided. This solver is used in Zingale & Katz (2015): From e05bc5c0f2ee193bd7847bac7e56b105e94fef93 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 10:19:39 -0400 Subject: [PATCH 21/25] add heading --- Util/exact_riemann/README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Util/exact_riemann/README.md b/Util/exact_riemann/README.md index 77dbd4d069..16fca8674b 100644 --- a/Util/exact_riemann/README.md +++ b/Util/exact_riemann/README.md @@ -1,3 +1,5 @@ +# Exact Riemann solver + This is a driver for the exact Riemann solver for a general equation of state. The main implementation is in Source/hydro. From b2b959257348f8932452a36e17c5249fd29c4562 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 10:58:35 -0400 Subject: [PATCH 22/25] add enum --- Source/hydro/riemann.cpp | 5 +++-- Source/hydro/riemann_constants.H | 8 ++++++++ Source/hydro/riemann_solvers.H | 4 ++-- 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index dc6eb0df98..ff9e876829 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -71,7 +71,8 @@ Castro::cmpflx_plus_godunov(const Box& bx, { - if (riemann_solver == 0 || riemann_solver == 1) { + if (riemann_solver == riemann_constants::TWO_SHOCK_CGF || + riemann_solver == riemann_constants::TWO_SHOCK_CG) { // approximate state Riemann solvers // first find the interface state on the current interface @@ -120,7 +121,7 @@ Castro::cmpflx_plus_godunov(const Box& bx, } } - } else if (riemann_solver == 2) { + } else if (riemann_solver == riemann_constants::HLLC) { // HLLC HLL::HLLC(i, j, k, idir, qm, qp, diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H index b29b92e2a5..eec8a6e501 100644 --- a/Source/hydro/riemann_constants.H +++ b/Source/hydro/riemann_constants.H @@ -6,6 +6,14 @@ using namespace amrex::literals; namespace riemann_constants { + + enum RiemannSolver : std::uint8_t { + TWO_SHOCK_CGF=0, + TWO_SHOCK_CG, + HLLC, + EXACT + }; + constexpr amrex::Real smlp1 = 1.e-10_rt; constexpr amrex::Real small = 1.e-8_rt; constexpr amrex::Real smallu = 1.e-12_rt; diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index 19b93287a7..76ccc0bc42 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -318,12 +318,12 @@ riemann_state(const int i, const int j, const int k, const int idir, // Solve Riemann problem - if (riemann_solver == 0) { + if (riemann_solver == riemann_constants::TWO_SHOCK_CGF) { // Colella, Glaz, & Ferguson solver TwoShock::riemannus(ql, qr, raux, qint); - } else if (riemann_solver == 1) { + } else if (riemann_solver == riemann_constants::TWO_SHOCK_CG) { // Colella & Glaz solver #ifndef RADIATION From e872e29f558801deadb5dcea6aa5f50caf28b608 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 15:46:09 -0400 Subject: [PATCH 23/25] a bit more work --- Source/hydro/exact_riemann_sample.H | 7 +------ Source/hydro/riemann.cpp | 3 ++- Source/hydro/riemann_solvers.H | 21 +++++++++++++++++++++ Util/exact_riemann/exact_riemann.H | 7 ++++++- 4 files changed, 30 insertions(+), 8 deletions(-) diff --git a/Source/hydro/exact_riemann_sample.H b/Source/hydro/exact_riemann_sample.H index c255ea1cce..a147763f6c 100644 --- a/Source/hydro/exact_riemann_sample.H +++ b/Source/hydro/exact_riemann_sample.H @@ -14,7 +14,7 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r, const amrex::Real ustar, const amrex::Real pstar, const amrex::Real W_l, const amrex::Real W_r, - const amrex::Real x, const amrex::Real xjump, const amrex::Real time, + const amrex::Real xi, amrex::Real& rho, amrex::Real& u, amrex::Real& p, amrex::Real* xn) { // get the initial sound speeds @@ -47,11 +47,6 @@ riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real // This follows from the discussion around C&G Eq. 15 - // compute xi = x/t -- this is the similarity variable for the - // solution - - amrex::Real xi = (x - xjump) / time; - // check which side of the contact we need to worry about amrex::Real chi = std::copysign(1.0_rt, xi - ustar); diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index ff9e876829..fe7efdc2fe 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -72,7 +72,8 @@ Castro::cmpflx_plus_godunov(const Box& bx, if (riemann_solver == riemann_constants::TWO_SHOCK_CGF || - riemann_solver == riemann_constants::TWO_SHOCK_CG) { + riemann_solver == riemann_constants::TWO_SHOCK_CG || + riemann_sovler == riemann_constants::EXACT) { // approximate state Riemann solvers // first find the interface state on the current interface diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index 76ccc0bc42..053947fbe6 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -330,6 +330,27 @@ riemann_state(const int i, const int j, const int k, const int idir, TwoShock::riemanncg(ql, qr, raux, qint); #endif + } else if (riemann_solver == riemann_constants::EXACT) { + // Exact Riemann solver + + // solve for the star state + + Real ustar{}; + Real pstar{}; + Real W_l{}; + Real W_r{}; + + auto niter = riemann_star_state(ql.rho, ql.un, ql.p, xn_l, + qr.rho, qr.un, qr.p, xn_r, + ustar, pstar, W_l, W_r); + + // sample the solution at xi = 0 + + + // we need to manually do the transverse sampling based on the + // contact + + #ifndef AMREX_USE_GPU } else { amrex::Error("ERROR: invalid value of riemann_solver"); diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index b1cf14e832..93e9140e52 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -78,10 +78,15 @@ exact_riemann() { amrex::Real rho, u, p, xn_s[NumSpec]; + // compute xi = x/t -- this is the similarity variable for the + // solution + + amrex::Real xi = (x - problem::xjump) / problem::t; + riemann_sample(problem::rho_l, problem::u_l, problem::p_l, xn, problem::rho_r, problem::u_r, problem::p_r, xn, ustar, pstar, W_l, W_r, - x, problem::xjump, problem::t, + xi, rho, u, p, xn_s); // get the thermodynamics for this state for output From 0fa676ce9ba5a33f93537f96d72806d4d4663248 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 17:59:04 -0400 Subject: [PATCH 24/25] more implementation --- Source/hydro/riemann.cpp | 2 +- Source/hydro/riemann_solvers.H | 58 +++++++++++++++++++++++++++++++--- 2 files changed, 55 insertions(+), 5 deletions(-) diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index fe7efdc2fe..2c2f284b80 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -73,7 +73,7 @@ Castro::cmpflx_plus_godunov(const Box& bx, if (riemann_solver == riemann_constants::TWO_SHOCK_CGF || riemann_solver == riemann_constants::TWO_SHOCK_CG || - riemann_sovler == riemann_constants::EXACT) { + riemann_solver == riemann_constants::EXACT) { // approximate state Riemann solvers // first find the interface state on the current interface diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index 053947fbe6..2b7f0e71fd 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -4,6 +4,8 @@ #include #include #include +#include +#include #ifdef HYBRID_MOMENTUM #include #endif @@ -335,10 +337,19 @@ riemann_state(const int i, const int j, const int k, const int idir, // solve for the star state - Real ustar{}; - Real pstar{}; - Real W_l{}; - Real W_r{}; + amrex::Real ustar{}; + amrex::Real pstar{}; + amrex::Real W_l{}; + amrex::Real W_r{}; + + amrex::Real xn_l[NumSpec]; + for (int n = 0; n < NumSpec; ++n) { + xn_l[n] = qm(i,j,k,QFS+n); + } + amrex::Real xn_r[NumSpec]; + for (int n = 0; n < NumSpec; ++n) { + xn_r[n] = qp(i,j,k,QFS+n); + } auto niter = riemann_star_state(ql.rho, ql.un, ql.p, xn_l, qr.rho, qr.un, qr.p, xn_r, @@ -346,10 +357,49 @@ riemann_state(const int i, const int j, const int k, const int idir, // sample the solution at xi = 0 + amrex::Real xi{}; + + amrex::Real rho{}; + amrex::Real u{}; + amrex::Real p{}; + amrex::Real xn_s[NumSpec]{}; + + riemann_sample(ql.rho, ql.un, ql.p, xn_l, + qr.rho, qr.un, qr.p, xn_r, + ustar, pstar, W_l, W_r, + xi, + rho, u, p, xn_s); + + // pack the RiemannState interface state + + qint.rho = rho; + qint.un = u; + qint.p = p; + + // get the energy + eos_rep_t eos_state; + eos_state.rho = qint.rho; + eos_state.T = castro::T_guess; + eos_state.p = qint.p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn_s[n]; + } + + eos(eos_input_rp, eos_state); + + qint.rhoe = qint.rho * eos_state.e; + qint.gamc = eos_state.gam1; // we need to manually do the transverse sampling based on the // contact + if (ustar > 0) { + qint.ut = ql.ut; + qint.utt = ql.utt; + } else { + qint.ut = qr.ut; + qint.utt = qr.utt; + } #ifndef AMREX_USE_GPU } else { From d060b0650f7fe2054824a4194c0f732832690ae9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jul 2024 18:20:47 -0400 Subject: [PATCH 25/25] fix compiler warning --- Source/hydro/riemann_solvers.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index 2b7f0e71fd..b547dc11ef 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -351,9 +351,9 @@ riemann_state(const int i, const int j, const int k, const int idir, xn_r[n] = qp(i,j,k,QFS+n); } - auto niter = riemann_star_state(ql.rho, ql.un, ql.p, xn_l, - qr.rho, qr.un, qr.p, xn_r, - ustar, pstar, W_l, W_r); + riemann_star_state(ql.rho, ql.un, ql.p, xn_l, + qr.rho, qr.un, qr.p, xn_r, + ustar, pstar, W_l, W_r); // sample the solution at xi = 0