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>

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

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<PulseProblem> {
 	static constexpr double mean_molecular_weight = mu;
@@ -234,6 +228,12 @@ auto problem_main() -> int
 	// Problem initialization
 	RadhydroSimulation<PulseProblem> 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)";

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<PulseProblem> {
 	static constexpr double mean_molecular_weight = mu;
@@ -239,6 +229,12 @@ auto problem_main() -> int
 	// Problem initialization
 	RadhydroSimulation<PulseProblem> 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<double> Tgas(nx);
 	std::vector<double> Vgas(nx);
 	std::vector<double> rhogas(nx);
+	std::vector<double> 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<PulseProblem>::gasDensity_index)[i];
 		const auto v_t = values.at(RadSystem<PulseProblem>::x1GasMomentum_index)[i] / rho_t;
 		const auto Egas = values.at(RadSystem<PulseProblem>::gasInternalEnergy_index)[i];
+		const auto flux_t = values.at(RadSystem<PulseProblem>::x1RadFlux_index)[i];
 		rhogas.at(i) = rho_t;
 		Trad.at(i) = Trad_t;
 		Tgas.at(i) = quokka::EOS<PulseProblem>::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<double> Tgas2(nx);
 	std::vector<double> Vgas2(nx);
 	std::vector<double> rhogas2(nx);
+	std::vector<double> 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<PulseProblem>::gasDensity_index)[i];
 		const auto v_t = values2.at(RadSystem<PulseProblem>::x1GasMomentum_index)[i] / rho_t;
 		const auto Egas = values2.at(RadSystem<PulseProblem>::gasInternalEnergy_index)[i];
+		const auto flux_t = values2.at(RadSystem<PulseProblem>::x1RadFlux_index)[i];
 		xs2.at(i) = x - drift;
 		rhogas2.at(index_) = rho_t;
 		Trad2.at(index_) = Trad_t;
 		Tgas2.at(index_) = quokka::EOS<PulseProblem>::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<std::string, std::string> 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

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<PulseProblem>::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<std::string, std::string> Tgas_args;
 	std::map<std::string, std::string> Texact_args;
 	std::map<std::string, std::string> 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

tests/MarshakAsymptotic.in

@@ -1,3 +1,5 @@
+cfl= = 10.0
+
 # *****************************************************************
 # Problem size and geometry
 # *****************************************************************

tests/RadhydroPulse.in

@@ -1,3 +1,7 @@
+kappa0 = 100.0
+v0_adv = 1.0e6
+max_time = 4.8e-5
+
 # *****************************************************************
 # Problem size and geometry
 # *****************************************************************

tests/RadhydroPulseDyn.in

@@ -1,3 +1,7 @@
+kappa0 = 500.0
+v0_adv = 3.0e7
+max_time = 4.8e-6
+
 # *****************************************************************
 # Problem size and geometry
 # *****************************************************************