From 54519af15fc1e323045cc9c1017195a3b2e452ff Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Wed, 10 Apr 2024 17:08:28 +1000 Subject: [PATCH] Update test parameters to produce all the figures in the IMEX paper (#600) ### Description Under the reproducible research philosophy which I believe in, I would like the reader to be able to reproduce all the figures and data from a paper with little effort. Therefore, I updated the tests used in the IMEX paper, setting some of the physical parameters as runtime parameters, so that the reader can reproduce the exact result in the paper by simply changing the config.in files, while the default parameters are kept moderate so that the automatic tests can run in a reasonable time. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [x] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .../test_radiation_marshak_asymptotic.cpp | 4 +- src/RadhydroPulse/test_radhydro_pulse.cpp | 22 ++++---- .../test_radhydro_pulse_dyn.cpp | 53 ++++++++++++++----- .../test_radhydro_uniform_advecting.cpp | 53 +++++++++++++------ tests/MarshakAsymptotic.in | 2 + tests/RadhydroPulse.in | 4 ++ tests/RadhydroPulseDyn.in | 4 ++ 7 files changed, 99 insertions(+), 43 deletions(-) diff --git a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index a529210c9..3c643fd4a 100644 --- a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -218,9 +218,9 @@ auto problem_main() -> int // Problem parameters const int max_timesteps = 1e6; - const double CFL_number = 0.9; + const double CFL_number = 10.0; const double initial_dt = 5.0e-12; // s - const double max_dt = 5.0e-12; // s + const double max_dt = 5.0; // s const double max_time = 10.0e-9; // s // const int nx = 60; // [18 == matches resolution of McClarren & Lowrie (2008)] // const double Lx = 0.66; // cm diff --git a/src/RadhydroPulse/test_radhydro_pulse.cpp b/src/RadhydroPulse/test_radhydro_pulse.cpp index 75989d2e8..f5a81548f 100644 --- a/src/RadhydroPulse/test_radhydro_pulse.cpp +++ b/src/RadhydroPulse/test_radhydro_pulse.cpp @@ -26,17 +26,11 @@ constexpr double width = 24.0; // cm, width of the pulse constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-10; constexpr double mu = 2.33 * C::m_u; constexpr double k_B = C::k_B; -constexpr double v0_nonadv = 0.; // non-advecting pulse -// static diffusion: tau = 2e3, beta = 3e-5, beta tau = 6e-2 -constexpr double kappa0 = 100.; // cm^2 g^-1 -constexpr double v0_adv = 1.0e6; // advecting pulse -constexpr double max_time = 4.8e-5; // max_time = 2.0 * width / v1; - -// dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 -// constexpr double kappa0 = 1000.; // cm^2 g^-1 -// constexpr double v0_adv = 1.0e8; // advecting pulse -// constexpr double max_time = 1.2e-4; // max_time = 2.0 * width / v1; +// Default parameters: static diffusion, tau = 2e3, beta = 3e-5, beta tau = 6e-2 +AMREX_GPU_MANAGED double kappa0 = 100.; // NOLINT +AMREX_GPU_MANAGED double v0_adv = 1.0e6; // NOLINT +// AMREX_GPU_MANAGED double max_time = 4.8e-5; // max_time = 2.0 * width / v1; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; @@ -234,6 +228,12 @@ auto problem_main() -> int // Problem initialization RadhydroSimulation sim(BCs_cc); + double max_time = 4.8e-5; + amrex::ParmParse pp; // NOLINT + pp.query("kappa0", kappa0); + pp.query("v0_adv", v0_adv); + pp.query("max_time", max_time); + sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; sim.radiationCflNumber_ = CFL_number; @@ -366,7 +366,7 @@ auto problem_main() -> int Trad_args["linestyle"] = "-."; Tgas_args["label"] = "Tgas (non-advecting)"; Tgas_args["linestyle"] = "--"; - matplotlibcpp::ylim(0.95e7, 1.6e7); + matplotlibcpp::ylim(0.95e7, 2.0e7); matplotlibcpp::plot(xs, Trad, Trad_args); matplotlibcpp::plot(xs, Tgas, Tgas_args); Trad_args["label"] = "Trad (advecting)"; diff --git a/src/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp b/src/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp index fb32d43b2..4b5120611 100644 --- a/src/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -27,21 +27,11 @@ constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-10; constexpr double mu = 2.33 * C::m_u; constexpr double k_B = C::k_B; -// Static diffusion: tau = 2e3, beta = 3e-5, beta tau = 6e-2 -// constexpr double kappa0 = 100.; // cm^2 g^-1 -// constexpr double v0_adv = 1.0e6; // advecting pulse -// constexpr double max_time = 4.8e-5; // max_time = 2.0 * width / v1; - -// Dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 -// constexpr double kappa0 = 1000.; // cm^2 g^-1 -// constexpr double v0_adv = 1.0e8; // advecting pulse -// constexpr double max_time = 4.8e-4; - -// Dynamic diffusion: tau = 1e4, beta = 1e-3, beta tau = 10. +// Default parameters: dynamic diffusion, tau = 1e4, beta = 1e-3, beta tau = 10. // Width of the pulse = sqrt(c max_time / kappa0) = 85 if max_time = 2.4e-4 -constexpr double kappa0 = 500.; // cm^2 g^-1 -constexpr double v0_adv = 3.0e7; // advecting pulse -constexpr double max_time = 4.8e-6; +AMREX_GPU_MANAGED double kappa0 = 500.; // NOLINT +AMREX_GPU_MANAGED double v0_adv = 3.0e7; // NOLINT +// AMREX_GPU_MANAGED double max_time = 4.8e-6; template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; @@ -239,6 +229,12 @@ auto problem_main() -> int // Problem initialization RadhydroSimulation sim(BCs_cc); + double max_time = 4.8e-6; + amrex::ParmParse pp; // NOLINT + pp.query("kappa0", kappa0); + pp.query("v0_adv", v0_adv); + pp.query("max_time", max_time); + sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; sim.radiationCflNumber_ = CFL_number; @@ -262,6 +258,7 @@ auto problem_main() -> int std::vector Tgas(nx); std::vector Vgas(nx); std::vector rhogas(nx); + std::vector flux(nx); for (int i = 0; i < nx; ++i) { amrex::Real const x = position[i]; @@ -271,10 +268,12 @@ auto problem_main() -> int const auto rho_t = values.at(RadSystem::gasDensity_index)[i]; const auto v_t = values.at(RadSystem::x1GasMomentum_index)[i] / rho_t; const auto Egas = values.at(RadSystem::gasInternalEnergy_index)[i]; + const auto flux_t = values.at(RadSystem::x1RadFlux_index)[i]; rhogas.at(i) = rho_t; Trad.at(i) = Trad_t; Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); Vgas.at(i) = 1e-5 * v_t; + flux.at(i) = flux_t; } // END OF PROBLEM 1 @@ -313,6 +312,7 @@ auto problem_main() -> int std::vector Tgas2(nx); std::vector Vgas2(nx); std::vector rhogas2(nx); + std::vector flux2(nx); for (int i = 0; i < nx; ++i) { int index_ = 0; @@ -335,11 +335,13 @@ auto problem_main() -> int const auto rho_t = values2.at(RadSystem::gasDensity_index)[i]; const auto v_t = values2.at(RadSystem::x1GasMomentum_index)[i] / rho_t; const auto Egas = values2.at(RadSystem::gasInternalEnergy_index)[i]; + const auto flux_t = values2.at(RadSystem::x1RadFlux_index)[i]; xs2.at(i) = x - drift; rhogas2.at(index_) = rho_t; Trad2.at(index_) = Trad_t; Tgas2.at(index_) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); Vgas2.at(index_) = 1e-5 * (v_t - v0_adv); + flux2.at(index_) = flux_t; } // END OF PROBLEM 2 @@ -435,6 +437,29 @@ auto problem_main() -> int } file.close(); + // plot radiation flux profile of the non-advecting pulse + matplotlibcpp::clf(); + std::map flux_args; + flux_args["label"] = "radiation flux (non-advecting)"; + flux_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, flux, flux_args); + flux_args["label"] = "radiation flux (advecting)"; + matplotlibcpp::plot(xs2, flux2, flux_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("flux (erg cm^-2 s^-1)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_dyndiff_flux.pdf"); + + // Save xs, flux, xs2, flux2 to csv file + file.open("radhydro_pulse_dyndiff_flux.csv"); + file << "xs,flux,xs2,flux2\n"; + for (size_t i = 0; i < xs.size(); ++i) { + file << std::scientific << std::setprecision(12) << xs[i] << "," << flux[i] << "," << xs2[i] << "," << flux2[i] << "\n"; + } + file.close(); + #endif // Cleanup and exit diff --git a/src/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp b/src/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp index f37c3447c..97c369c69 100644 --- a/src/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp +++ b/src/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp @@ -12,24 +12,39 @@ struct PulseProblem { }; // dummy type to allow compile-type polymorphism via template specialization +constexpr double c = 1.0e8; +// model 0 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-4 * c; +// constexpr double kappa0 = 1.0e4; // dx = 1, tau = kappa0 * dx = 1e4 +// constexpr double chat = 1.0e7; +// model 1 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-4 * c; +// constexpr double kappa0 = 1.0e4; // dx = 1, tau = kappa0 * dx = 1e4 +// constexpr double chat = 1.0e8; +// model 2 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-2 * c; +// constexpr double kappa0 = 1.0e5; +// constexpr double chat = 1.0e8; +// model 3 constexpr int beta_order_ = 2; // order of beta in the radiation four-force +constexpr double v0 = 1e-2 * c; +constexpr double kappa0 = 1.0e5; +constexpr double chat = 1.0e8; constexpr double T0 = 1.0; // temperature constexpr double rho0 = 1.0; // matter density constexpr double a_rad = 1.0; -constexpr double c = 1.0; -constexpr double chat = c; constexpr double mu = 1.0; constexpr double k_B = 1.0; // static diffusion, beta = 1e-4, tau_cell = kappa0 * dx = 100, beta tau_cell = 1e-2 // constexpr double kappa0 = 100.; // cm^2 g^-1 // constexpr double v0 = 1.0e-4 * c; // advecting pulse -// constexpr double max_time = 1.0 / v0; // dynamic diffusion, beta = 1e-3, tau = kappa0 * dx = 1e5, beta tau = 100 -constexpr double kappa0 = 1.0e4; // dx = 1, tau = kappa0 * dx = 1e4 -constexpr double v0 = 1e-3 * c; // beta = 1e-3 constexpr double max_time = 10.0 / v0; constexpr double Erad0 = a_rad * T0 * T0 * T0 * T0; @@ -86,7 +101,10 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(qu double erad = NAN; double frad = NAN; - if constexpr (beta_order_ <= 1) { + if constexpr (beta_order_ == 0) { + erad = Erad0; + frad = 0.0; + } else if constexpr (beta_order_ == 1) { erad = Erad0; frad = 4. / 3. * v0 * Erad0; } else if constexpr (beta_order_ == 2) { @@ -125,8 +143,9 @@ auto problem_main() -> int // of order 10^5. // Problem parameters - const int max_timesteps = 1e5; - const double CFL_number = 0.8; + const int max_timesteps = 1e6; + const double CFL_number_gas = 0.8; + const double CFL_number_rad = 8.0; const double max_dt = 1.0; @@ -145,8 +164,8 @@ auto problem_main() -> int sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; - sim.radiationCflNumber_ = CFL_number; - sim.cflNumber_ = CFL_number; + sim.radiationCflNumber_ = CFL_number_rad; + sim.cflNumber_ = CFL_number_gas; sim.maxDt_ = max_dt; sim.maxTimesteps_ = max_timesteps; sim.plotfileInterval_ = -1; @@ -218,22 +237,25 @@ auto problem_main() -> int std::map Tgas_args; std::map Texact_args; std::map Tradexact_args; - Trad_args["label"] = "radiation temperature"; + Trad_args["label"] = "radiation (numerical)"; Trad_args["linestyle"] = "-"; - Tradexact_args["label"] = "radiation temperature (exact)"; + Tradexact_args["label"] = "radiation (exact)"; Tradexact_args["linestyle"] = "--"; - Tgas_args["label"] = "gas temperature"; + Tgas_args["label"] = "gas (numerical)"; Tgas_args["linestyle"] = "-"; - Texact_args["label"] = "gas temperature (exact)"; + Texact_args["label"] = "gas (exact)"; Texact_args["linestyle"] = "--"; matplotlibcpp::plot(xs, Trad, Trad_args); matplotlibcpp::plot(xs, Trad_exact, Tradexact_args); matplotlibcpp::plot(xs, Tgas, Tgas_args); matplotlibcpp::plot(xs, Tgas_exact, Texact_args); - matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::xlabel("x (dimensionless)"); matplotlibcpp::ylabel("temperature (dimensionless)"); matplotlibcpp::legend(); matplotlibcpp::title(fmt::format("time ct = {:.4g}", sim.tNew_[0] * c)); + if constexpr (beta_order_ == 1) { + matplotlibcpp::ylim(1.0 - 1.0e-7, 1.0 + 1.0e-7); + } matplotlibcpp::tight_layout(); matplotlibcpp::save("./radhydro_uniform_advecting_temperature_dimensionless.pdf"); @@ -252,7 +274,6 @@ auto problem_main() -> int matplotlibcpp::title(fmt::format("time ct = {:.4g}", sim.tNew_[0] * c)); matplotlibcpp::tight_layout(); matplotlibcpp::save("./radhydro_uniform_advecting_velocity_dimensionless.pdf"); - #endif // Cleanup and exit diff --git a/tests/MarshakAsymptotic.in b/tests/MarshakAsymptotic.in index a87383dec..4bd767cb1 100644 --- a/tests/MarshakAsymptotic.in +++ b/tests/MarshakAsymptotic.in @@ -1,3 +1,5 @@ +cfl= = 10.0 + # ***************************************************************** # Problem size and geometry # ***************************************************************** diff --git a/tests/RadhydroPulse.in b/tests/RadhydroPulse.in index 51ba0ecfb..154f44d9f 100644 --- a/tests/RadhydroPulse.in +++ b/tests/RadhydroPulse.in @@ -1,3 +1,7 @@ +kappa0 = 100.0 +v0_adv = 1.0e6 +max_time = 4.8e-5 + # ***************************************************************** # Problem size and geometry # ***************************************************************** diff --git a/tests/RadhydroPulseDyn.in b/tests/RadhydroPulseDyn.in index 3931464c2..48daab6dc 100644 --- a/tests/RadhydroPulseDyn.in +++ b/tests/RadhydroPulseDyn.in @@ -1,3 +1,7 @@ +kappa0 = 500.0 +v0_adv = 3.0e7 +max_time = 4.8e-6 + # ***************************************************************** # Problem size and geometry # *****************************************************************