From 7a597a57579780a8e0e25ee0220296de479ff4ea Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Wed, 29 May 2024 03:23:51 +1000 Subject: [PATCH] Fix radiation transport error in Y and Z directions. (#633) ### Description Fix bug #632 . This is a dumb mistake I made on using `AMREX_D_TERM` in `AddFluxesRK2`, causing a wrong calculation of the HLL flux in 2D/3D cases. The mistake was introduced in the "Use IMEX PD-ARS ..." commit (78c5970857430e155885555bdfaa30a6ef4946e4) 6 months ago. Because of this mistake, all 2D/3D simulations since then with non-zero y/z-component radiation fluxes are wrong. To avoid mistakes like this in the future, I added a new test of radiation streaming in the Y direction. I plan to add a better test (say, the RadShock test) in both Y and Z direction, with the following grid configuration: [4, 128, 4] and [4, 4, 128]. Or, alternatively, we should have a non-trivial 3D test with an exact solution. ### Related issues Fixed #632 ### 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/CMakeLists.txt | 1 + src/RadStreamingY/CMakeLists.txt | 9 + .../test_radiation_streaming_y.cpp | 270 ++++++++++++++++++ .../test_radiation_streaming_y.hpp | 24 ++ src/radiation_system.hpp | 4 +- tests/RadStreamingY.in | 24 ++ 6 files changed, 330 insertions(+), 2 deletions(-) create mode 100644 src/RadStreamingY/CMakeLists.txt create mode 100644 src/RadStreamingY/test_radiation_streaming_y.cpp create mode 100644 src/RadStreamingY/test_radiation_streaming_y.hpp create mode 100644 tests/RadStreamingY.in diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 914bdfa6e..5ea736290 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -164,6 +164,7 @@ add_subdirectory(RadMatterCouplingRSLA) add_subdirectory(RadPulse) add_subdirectory(RadShadow) add_subdirectory(RadStreaming) +add_subdirectory(RadStreamingY) add_subdirectory(RadSuOlson) add_subdirectory(RadTophat) add_subdirectory(RadTube) diff --git a/src/RadStreamingY/CMakeLists.txt b/src/RadStreamingY/CMakeLists.txt new file mode 100644 index 000000000..2ddc62d1e --- /dev/null +++ b/src/RadStreamingY/CMakeLists.txt @@ -0,0 +1,9 @@ +if (AMReX_SPACEDIM GREATER_EQUAL 2) + add_executable(test_radiation_streaming_y test_radiation_streaming_y.cpp ../fextract.cpp ${QuokkaObjSources}) + + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radiation_streaming_y) + endif(AMReX_GPU_BACKEND MATCHES "CUDA") + + add_test(NAME RadiationStreamingY COMMAND test_radiation_streaming_y RadStreamingY.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/RadStreamingY/test_radiation_streaming_y.cpp b/src/RadStreamingY/test_radiation_streaming_y.cpp new file mode 100644 index 000000000..7f2c08505 --- /dev/null +++ b/src/RadStreamingY/test_radiation_streaming_y.cpp @@ -0,0 +1,270 @@ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radiation_streaming.cpp +/// \brief Defines a test problem for radiation in the free-streaming regime. +/// + +#include "test_radiation_streaming_y.hpp" +#include "AMReX.H" +#include "AMReX_BLassert.H" +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "valarray.hpp" + +struct StreamingProblem { +}; + +constexpr int direction = 1; + +constexpr double initial_Erad = 1.0e-5; +constexpr double initial_Egas = 1.0e-5; +constexpr double c = 1.0; // speed of light +constexpr double chat = c; // reduced speed of light +constexpr double kappa0 = 1.0e-10; // opacity +constexpr double rho = 1.0; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = 1.0; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = false; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = 1; // number of radiation groups +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = 1.0; + static constexpr double Erad_floor = initial_Erad; + static constexpr int beta_order = 0; +}; + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray +{ + quokka::valarray kappaPVec{}; + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = kappa0; + } + return kappaPVec; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) + -> quokka::valarray +{ + return ComputePlanckOpacity(0.0, 0.0); +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const auto Erad0 = initial_Erad; + const auto Egas0 = initial_Egas; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + state_cc(i, j, k, RadSystem::radEnergy_index) = Erad0; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::x2RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas0; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, int /*dcomp*/, + int /*numcomp*/, amrex::GeometryData const &geom, const amrex::Real /*time*/, + const amrex::BCRec * /*bcr*/, int /*bcomp*/, int /*orig_comp*/) +{ +#if (AMREX_SPACEDIM == 1) + auto i = iv.toArray()[0]; + int j = 0; + int k = 0; +#endif +#if (AMREX_SPACEDIM == 2) + auto [i, j] = iv.toArray(); + int k = 0; +#endif +#if (AMREX_SPACEDIM == 3) + auto [i, j, k] = iv.toArray(); +#endif + + amrex::Box const &box = geom.Domain(); + amrex::GpuArray lo = box.loVect3d(); + amrex::GpuArray hi = box.hiVect3d(); + + if constexpr (direction == 0) { + if (i < lo[0]) { + // streaming inflow boundary + const double Erad = 1.0; + const double Frad = c * Erad; + + // x1 left side boundary (Marshak) + consVar(i, j, k, RadSystem::radEnergy_index) = Erad; + consVar(i, j, k, RadSystem::x1RadFlux_index) = Frad; + consVar(i, j, k, RadSystem::x2RadFlux_index) = 0.; + consVar(i, j, k, RadSystem::x3RadFlux_index) = 0.; + } else if (i >= hi[0]) { + // right-side boundary -- constant + const double Erad = initial_Erad; + consVar(i, j, k, RadSystem::radEnergy_index) = Erad; + consVar(i, j, k, RadSystem::x1RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x2RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index) = 0; + } + } else { + if (j < lo[1]) { + // streaming inflow boundary + const double Erad = 1.0; + const double Frad = c * Erad; + + // x1 left side boundary (Marshak) + consVar(i, j, k, RadSystem::radEnergy_index) = Erad; + consVar(i, j, k, RadSystem::x1RadFlux_index) = 0.; + consVar(i, j, k, RadSystem::x2RadFlux_index) = Frad; + consVar(i, j, k, RadSystem::x3RadFlux_index) = 0.; + } else if (j >= hi[1]) { + // right-side boundary -- constant + const double Erad = initial_Erad; + consVar(i, j, k, RadSystem::radEnergy_index) = Erad; + consVar(i, j, k, RadSystem::x1RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x2RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index) = 0; + } + } + + // gas boundary conditions are the same everywhere + const double Egas = initial_Egas; + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas; + consVar(i, j, k, RadSystem::gasDensity_index) = rho; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; +} + +auto problem_main() -> int +{ + // Problem parameters + // const int nx = 1000; + // const double Lx = 1.0; + const double CFL_number = 0.8; + const double dt_max = 1e-2; + double tmax = 1.0; + const int max_timesteps = 5000; + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + // assert at compile time + if constexpr (direction == 0) { + BCs_cc[n].setLo(0, amrex::BCType::ext_dir); // Dirichlet x1 + BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1 + BCs_cc[n].setLo(1, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(1, amrex::BCType::int_dir); + } else { + BCs_cc[n].setLo(0, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(0, amrex::BCType::int_dir); + BCs_cc[n].setLo(1, amrex::BCType::ext_dir); // Dirichlet x1 + BCs_cc[n].setHi(1, amrex::BCType::foextrap); // extrapolate x1 + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + // read tmax from inputs file + amrex::ParmParse pp; // NOLINT + pp.query("max_time", tmax); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = tmax; + sim.radiationCflNumber_ = CFL_number; + sim.maxDt_ = dt_max; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), direction, 0.0); + const int nx = static_cast(position.size()); + + // compute error norm + std::vector erad(nx); + std::vector erad_exact(nx); + std::vector xs(nx); + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + erad_exact.at(i) = (x <= chat * tmax) ? 1.0 : 0.0; + erad.at(i) = values.at(RadSystem::radEnergy_index)[i]; + } + + double err_norm = 0.; + double sol_norm = 0.; + for (int i = 0; i < nx; ++i) { + err_norm += std::abs(erad[i] - erad_exact[i]); + sol_norm += std::abs(erad_exact[i]); + } + + const double rel_err_norm = err_norm / sol_norm; + const double rel_err_tol = 0.05; + int status = 1; + if (rel_err_norm < rel_err_tol) { + status = 0; + } + amrex::Print() << "Relative L1 norm = " << rel_err_norm << std::endl; + +#ifdef HAVE_PYTHON + // Plot results + matplotlibcpp::clf(); + matplotlibcpp::ylim(-0.1, 1.1); + + std::map erad_args; + std::map erad_exact_args; + erad_args["label"] = "numerical solution"; + erad_exact_args["label"] = "exact solution"; + erad_exact_args["linestyle"] = "--"; + matplotlibcpp::plot(xs, erad, erad_args); + matplotlibcpp::plot(xs, erad_exact, erad_exact_args); + + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("t = {:.4f}", sim.tNew_[0])); + if constexpr (direction == 0) { + matplotlibcpp::save("./radiation_streaming_x.pdf"); + } else { + matplotlibcpp::save("./radiation_streaming_y.pdf"); + } +#endif // HAVE_PYTHON + + // Cleanup and exit + amrex::Print() << "Finished." << std::endl; + return status; +} diff --git a/src/RadStreamingY/test_radiation_streaming_y.hpp b/src/RadStreamingY/test_radiation_streaming_y.hpp new file mode 100644 index 000000000..11f54b077 --- /dev/null +++ b/src/RadStreamingY/test_radiation_streaming_y.hpp @@ -0,0 +1,24 @@ +#ifndef TEST_RADIATION_STREAMING_Y_HPP_ // NOLINT +#define TEST_RADIATION_STREAMING_Y_HPP_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radiation_streaming.hpp +/// \brief Defines a test problem for radiation in the free-streaming regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include + +// internal headers + +#include "radiation_system.hpp" + +// function definitions + +#endif // TEST_RADIATION_STREAMING_Y_HPP_ diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index b3fb6e27b..4292b72c2 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -508,8 +508,8 @@ void RadSystem::AddFluxesRK2(array_t &U_new, arrayconst_t &U0, arrayc const double FzU_1 = (dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n)); #endif // save results in cons_new - cons_new[n] = (1.0 - IMEX_a32) * U_0 + IMEX_a32 * U_1 + ((0.5 - IMEX_a32) * AMREX_D_TERM(FxU_0, +FyU_0, +FzU_0)) + - (0.5 * AMREX_D_TERM(FxU_1, +FyU_1, +FzU_1)); + cons_new[n] = (1.0 - IMEX_a32) * U_0 + IMEX_a32 * U_1 + ((0.5 - IMEX_a32) * (AMREX_D_TERM(FxU_0, +FyU_0, +FzU_0))) + + (0.5 * (AMREX_D_TERM(FxU_1, +FyU_1, +FzU_1))); } if (!isStateValid(cons_new)) { diff --git a/tests/RadStreamingY.in b/tests/RadStreamingY.in new file mode 100644 index 000000000..c30fcd360 --- /dev/null +++ b/tests/RadStreamingY.in @@ -0,0 +1,24 @@ +max_time = 0.2 + +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 1 0 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 4 100 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 0 \ No newline at end of file