From 3455fc8ef2ce93033df9d7be11651e427a1f49a1 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 3 Feb 2024 07:54:24 +1100 Subject: [PATCH] Multigroup advecting pulse test (#486) ### Description Define the multigroup, variable opacity advecting pulse test. See Sec. 4.1.6 of the CASTRO III paper (Zhang et al. 2013). The code perfectly reproduces Figure 6 of the paper. We use the original Newton-Raphson method on the (Egas, Erad) basis. ### Related issues N/A ### 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). - [ ] 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> --- .ci/azure-pipelines-aarch64-debug.yml | 2 +- src/CMakeLists.txt | 1 + src/RadhydroPulse/CMakeLists.txt | 2 +- src/RadhydroPulseMG/CMakeLists.txt | 9 + .../test_radhydro_pulse_MG.cpp | 523 ++++++++++++++++++ .../test_radhydro_pulse_MG.hpp | 21 + .../test_radhydro_shock_multigroup.cpp | 3 +- src/planck_integral.hpp | 3 + src/radiation_system.hpp | 172 +++--- 9 files changed, 662 insertions(+), 74 deletions(-) create mode 100644 src/RadhydroPulseMG/CMakeLists.txt create mode 100644 src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp create mode 100644 src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp diff --git a/.ci/azure-pipelines-aarch64-debug.yml b/.ci/azure-pipelines-aarch64-debug.yml index 3d53e9d87..6da3152a7 100644 --- a/.ci/azure-pipelines-aarch64-debug.yml +++ b/.ci/azure-pipelines-aarch64-debug.yml @@ -30,7 +30,7 @@ jobs: - task: CMake@1 displayName: 'Run CTest' inputs: - cmakeArgs: '-E chdir . ctest -j 2 -T Test --output-on-failure' + cmakeArgs: '-E chdir . ctest -j 2 -T Test -E "RadhydroPulseMG*" --output-on-failure' - task: PublishTestResults@2 inputs: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e38854230..689728792 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -172,6 +172,7 @@ add_subdirectory(RadhydroShockCGS) add_subdirectory(RadhydroShockMultigroup) add_subdirectory(RadhydroUniformAdvecting) add_subdirectory(RadhydroPulse) +add_subdirectory(RadhydroPulseMG) add_subdirectory(ODEIntegration) add_subdirectory(Cooling) diff --git a/src/RadhydroPulse/CMakeLists.txt b/src/RadhydroPulse/CMakeLists.txt index c77a58d81..0cc17b65f 100644 --- a/src/RadhydroPulse/CMakeLists.txt +++ b/src/RadhydroPulse/CMakeLists.txt @@ -6,4 +6,4 @@ if (AMReX_SPACEDIM EQUAL 1) endif(AMReX_GPU_BACKEND MATCHES "CUDA") add_test(NAME RadhydroPulse COMMAND test_radhydro_pulse RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) -endif() \ No newline at end of file +endif() diff --git a/src/RadhydroPulseMG/CMakeLists.txt b/src/RadhydroPulseMG/CMakeLists.txt new file mode 100644 index 000000000..e62068b4f --- /dev/null +++ b/src/RadhydroPulseMG/CMakeLists.txt @@ -0,0 +1,9 @@ +if (AMReX_SPACEDIM EQUAL 1) + add_executable(test_radhydro_pulse_MG test_radhydro_pulse_MG.cpp ../fextract.cpp ${QuokkaObjSources}) + + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_pulse_MG) + endif(AMReX_GPU_BACKEND MATCHES "CUDA") + + add_test(NAME RadhydroPulseMG COMMAND test_radhydro_pulse_MG RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp new file mode 100644 index 000000000..23ca8a5c0 --- /dev/null +++ b/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp @@ -0,0 +1,523 @@ +/// \file test_radhydro_pulse_MG.cpp +/// \brief Defines a test problem for multigroup radiation in the diffusion regime with advection by gas. +/// + +#include "test_radhydro_pulse_MG.hpp" +#include "AMReX_BC_TYPES.H" +#include "AMReX_Print.H" +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "physics_info.hpp" + +struct PulseProblem { +}; // dummy type to allow compile-type polymorphism via template specialization +struct AdvPulseProblem { +}; + +constexpr int isconst = 0; +// constexpr int n_groups_ = 1; +// constexpr amrex::GpuArray rad_boundaries_{0., inf}; +constexpr int n_groups_ = 2; +constexpr amrex::GpuArray rad_boundaries_{1e15, 1e17, 1e19}; +// constexpr int n_groups_ = 8; +// constexpr amrex::GpuArray rad_boundaries_{1e15, 3.16e15, 1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19}; +// constexpr int n_groups_ = 64; +// constexpr amrex::GpuArray rad_boundaries_{1.00000000e+15, 1.15478198e+15, 1.33352143e+15, 1.53992653e+15, +// 1.77827941e+15, 2.05352503e+15, 2.37137371e+15, 2.73841963e+15, +// 3.16227766e+15, 3.65174127e+15, 4.21696503e+15, 4.86967525e+15, +// 5.62341325e+15, 6.49381632e+15, 7.49894209e+15, 8.65964323e+15, +// 1.00000000e+16, 1.15478198e+16, 1.33352143e+16, 1.53992653e+16, +// 1.77827941e+16, 2.05352503e+16, 2.37137371e+16, 2.73841963e+16, +// 3.16227766e+16, 3.65174127e+16, 4.21696503e+16, 4.86967525e+16, +// 5.62341325e+16, 6.49381632e+16, 7.49894209e+16, 8.65964323e+16, +// 1.00000000e+17, 1.15478198e+17, 1.33352143e+17, 1.53992653e+17, +// 1.77827941e+17, 2.05352503e+17, 2.37137371e+17, 2.73841963e+17, +// 3.16227766e+17, 3.65174127e+17, 4.21696503e+17, 4.86967525e+17, +// 5.62341325e+17, 6.49381632e+17, 7.49894209e+17, 8.65964323e+17, +// 1.00000000e+18, 1.15478198e+18, 1.33352143e+18, 1.53992653e+18, +// 1.77827941e+18, 2.05352503e+18, 2.37137371e+18, 2.73841963e+18, +// 3.16227766e+18, 3.65174127e+18, 4.21696503e+18, 4.86967525e+18, +// 5.62341325e+18, 6.49381632e+18, 7.49894209e+18, 8.65964323e+18, +// 1.00000000e+19}; + +constexpr double kappa0 = 180.; // cm^2 g^-1 +constexpr double kappa_const = 100.; // cm^2 g^-1 +constexpr double T0 = 1.0e7; // K (temperature) +constexpr double T1 = 2.0e7; // K (temperature) +constexpr double rho0 = 1.2; // g cm^-3 (matter density) +constexpr double a_rad = C::a_rad; +constexpr double c = C::c_light; // speed of light (cgs) +constexpr double chat = c; +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 h_planck = C::hplanck; +constexpr double k_B = C::k_B; +constexpr double v0_nonadv = 0.; // non-advecting pulse + +// static diffusion: (for single group) tau = 2e3, beta = 3e-5, beta tau = 6e-2 +constexpr double v0_adv = 1.0e6; // advecting pulse +constexpr double max_time = 2.4e-5; // max_time = 1.0 * width / v1; +// 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; + +constexpr double T_ref = T0; +constexpr double nu_ref = 1.0e18; // Hz +constexpr double coeff_ = h_planck * nu_ref / (k_B * T_ref); + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + 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 = n_groups_; +}; +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + 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 = n_groups_; +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr bool compute_v_over_c_terms = true; + static constexpr double energy_unit = h_planck; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; +}; +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr bool compute_v_over_c_terms = true; + static constexpr double energy_unit = h_planck; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; +}; + +AMREX_GPU_HOST_DEVICE +auto compute_initial_Tgas(const double x) -> double +{ + // compute temperature profile for Gaussian radiation pulse + const double sigma = width; + return T0 + (T1 - T0) * std::exp(-x * x / (2.0 * sigma * sigma)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_exact_rho(const double x) -> double +{ + // compute density profile for Gaussian radiation pulse + auto T = compute_initial_Tgas(x); + return rho0 * T0 / T + (a_rad * mu / 3. / k_B) * (std::pow(T0, 4) / T - std::pow(T, 3)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_kappa(const double nu, const double Tgas) -> double +{ + // cm^-1 + auto T_ = Tgas / T_ref; + auto nu_ = nu / nu_ref; + return kappa0 * std::pow(T_, -0.5) * std::pow(nu_, -3) * (1.0 - std::exp(-coeff_ * nu_ / T_)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_repres_nu() -> quokka::valarray +{ + // return the geometrical mean as the representative frequency for each group + quokka::valarray nu_rep{}; + if constexpr (n_groups_ == 1) { + nu_rep[0] = nu_ref; + } else { + for (int g = 0; g < n_groups_; ++g) { + nu_rep[g] = std::sqrt(rad_boundaries_[g] * rad_boundaries_[g + 1]); + } + } + return nu_rep; +} + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + quokka::valarray kappaPVec{}; + auto nu_rep = compute_repres_nu(); + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = compute_kappa(nu_rep[g], Tgas) / rho; + } + if constexpr (isconst == 1) { + kappaPVec.fillin(kappa_const); + } + return kappaPVec; +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return RadSystem::ComputePlanckOpacity(rho, Tgas); +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) + -> quokka::valarray +{ + quokka::valarray opacity_deriv{}; + const auto nu_rep = compute_repres_nu(); + const auto T = Tgas / T0; + for (int g = 0; g < nGroups_; ++g) { + const auto nu = nu_rep[g] / nu_ref; + opacity_deriv[g] = + 1. / T0 * kappa0 * (-0.5 * std::pow(T, -1.5) * (1. - std::exp(-coeff_ * nu / T)) - coeff_ * std::pow(T, -2.5) * std::exp(-coeff_ * nu / T)); + opacity_deriv[g] /= rho; + } + if constexpr (isconst == 1) { + opacity_deriv.fillin(0.0); + } + return opacity_deriv; +} +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) + -> quokka::valarray +{ + return RadSystem::ComputePlanckOpacityTempDerivative(rho, Tgas); +} + +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double +{ + return (1. / 3.); // Eddington approximation +} +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double +{ + return (1. / 3.); // Eddington approximation +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + amrex::GpuArray prob_hi = grid_elem.prob_hi_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); + + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + const double Trad = compute_initial_Tgas(x - x0); + const double rho = compute_exact_rho(x - x0); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double v0 = v0_nonadv; + + auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + amrex::GpuArray prob_hi = grid_elem.prob_hi_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); + + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + const double Trad = compute_initial_Tgas(x - x0); + const double rho = compute_exact_rho(x - x0); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double v0 = v0_adv; + + auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +auto problem_main() -> int +{ + // This problem is a test of radiation diffusion plus advection by gas. + // This makes this problem a stringent test of the radiation advection + // in the diffusion limit. + + // Problem parameters + const int64_t max_timesteps = 1e8; + const double CFL_number = 0.8; + // const int nx = 32; + + const double max_dt = 1e-3; // t_cr = 2 cm / cs = 7e-8 s + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + // periodic boundary condition in the x-direction will not work + BCs_cc[n].setLo(0, amrex::BCType::foextrap); // extrapolate + BCs_cc[n].setHi(0, amrex::BCType::foextrap); + for (int i = 1; i < AMREX_SPACEDIM; ++i) { + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); + } + } + + // Problem 1: non-advecting pulse + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + sim.radiationCflNumber_ = CFL_number; + sim.cflNumber_ = CFL_number; + sim.maxDt_ = max_dt; + 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), 0, 0.0); + const int nx = static_cast(position.size()); + amrex::GpuArray prob_lo = sim.geom[0].ProbLoArray(); + amrex::GpuArray prob_hi = sim.geom[0].ProbHiArray(); + + std::vector xs(nx); + std::vector Trad(nx); + std::vector Tgas(nx); + std::vector Vgas(nx); + std::vector rhogas(nx); + + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + // const auto Erad_t = values.at(RadSystem::radEnergy_index)[i]; + double Erad_t = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + 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]; + 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; + } + // END OF PROBLEM 1 + + // Problem 2: advecting pulse + + // Problem initialization + RadhydroSimulation sim2(BCs_cc); + + sim2.radiationReconstructionOrder_ = 3; // PPM + sim2.stopTime_ = max_time; + sim2.radiationCflNumber_ = CFL_number; + sim2.maxDt_ = max_dt; + sim2.maxTimesteps_ = max_timesteps; + sim2.plotfileInterval_ = -1; + + // initialize + sim2.setInitialConditions(); + + // evolve + sim2.evolve(); + + // read output variables + auto [position2, values2] = fextract(sim2.state_new_cc_[0], sim2.Geom(0), 0, 0.0); + prob_lo = sim2.geom[0].ProbLoArray(); + prob_hi = sim2.geom[0].ProbHiArray(); + // compute the pixel size + const double dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx); + const double move = v0_adv * sim2.tNew_[0]; + const int n_p = static_cast(move / dx); + const int half = static_cast(nx / 2.0); + const double drift = move - static_cast(n_p) * dx; + const int shift = n_p - static_cast((n_p + half) / nx) * nx; + + std::vector xs2(nx); + std::vector Trad2(nx); + std::vector Tgas2(nx); + std::vector Vgas2(nx); + std::vector rhogas2(nx); + + for (int i = 0; i < nx; ++i) { + int index_ = 0; + if (shift >= 0) { + if (i < shift) { + index_ = nx - shift + i; + } else { + index_ = i - shift; + } + } else { + if (i <= nx - 1 + shift) { + index_ = i - shift; + } else { + index_ = i - (nx + shift); + } + } + const amrex::Real x = position2[i]; + // const auto Erad_t = values2.at(RadSystem::radEnergy_index)[i]; + double Erad_t = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values2.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + 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]; + 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); + } + // END OF PROBLEM 2 + + // compute error norm + double err_norm = 0.; + double sol_norm = 0.; + for (size_t i = 0; i < xs2.size(); ++i) { + err_norm += std::abs(Tgas[i] - Trad[i]); + err_norm += std::abs(Trad2[i] - Trad[i]); + err_norm += std::abs(Tgas2[i] - Trad[i]); + sol_norm += std::abs(Trad[i]) * 3.0; + } + const double error_tol = 0.008; + const double rel_error = err_norm / sol_norm; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + +#ifdef HAVE_PYTHON + // plot temperature + matplotlibcpp::clf(); + std::map Trad_args; + std::map Tgas_args; + Trad_args["label"] = "Trad (non-advecting)"; + Trad_args["linestyle"] = "-."; + Tgas_args["label"] = "Tgas (non-advecting)"; + Tgas_args["linestyle"] = "--"; + matplotlibcpp::plot(xs, Trad, Trad_args); + matplotlibcpp::plot(xs, Tgas, Tgas_args); + Trad_args["label"] = "Trad (advecting)"; + Tgas_args["label"] = "Tgas (advecting)"; + matplotlibcpp::plot(xs2, Trad2, Trad_args); + matplotlibcpp::plot(xs2, Tgas2, Tgas_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("temperature (K)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_temperature.pdf"); + + // plot gas density profile + matplotlibcpp::clf(); + std::map rho_args; + rho_args["label"] = "gas density (non-advecting)"; + rho_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, rhogas, rho_args); + rho_args["label"] = "gas density (advecting))"; + matplotlibcpp::plot(xs2, rhogas2, rho_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("density (g cm^-3)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_density.pdf"); + + // plot gas velocity profile + matplotlibcpp::clf(); + std::map vgas_args; + vgas_args["label"] = "gas velocity (non-advecting)"; + vgas_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, Vgas, vgas_args); + vgas_args["label"] = "gas velocity (advecting)"; + matplotlibcpp::plot(xs2, Vgas2, vgas_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("velocity (km s^-1)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_velocity.pdf"); + +#endif + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + return status; +} diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp new file mode 100644 index 000000000..eed311866 --- /dev/null +++ b/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp @@ -0,0 +1,21 @@ +#ifndef TEST_RADHYDRO_PULSE_MG_HPP_ // NOLINT +#define TEST_RADHYDRO_PULSE_MG_HPP_ +/// \file test_radiation_marshak.hpp +/// \brief Defines a test problem for radiation in the static diffusion regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// function definitions + +#endif // TEST_RADHYDRO_PULSE_MG_HPP_ diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index 1640fb193..1a21b2dd1 100644 --- a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -40,6 +40,7 @@ constexpr double v1 = 1.73e7; // cm s^-1 [1.72875e7] constexpr double chat = 10.0 * (v0 + c_s0); // reduced speed of light constexpr double Erad0 = a_rad * (T0 * T0 * T0 * T0); // erg cm^-3 +constexpr double Erad_floor_ = Erad0 * 1e-12; constexpr double Egas0 = rho0 * c_v * T0; // erg cm^-3 constexpr double Erad1 = a_rad * (T1 * T1 * T1 * T1); // erg cm^-3 constexpr double Egas1 = rho1 * c_v * T1; // erg cm^-3 @@ -63,7 +64,7 @@ template <> struct RadSystem_Traits { static constexpr double c_light = c; static constexpr double c_hat = chat; static constexpr double radiation_constant = a_rad; - static constexpr double Erad_floor = 0.; + static constexpr double Erad_floor = Erad_floor_; static constexpr bool compute_v_over_c_terms = true; static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{ diff --git a/src/planck_integral.hpp b/src/planck_integral.hpp index deb22e8ac..56d53a337 100644 --- a/src/planck_integral.hpp +++ b/src/planck_integral.hpp @@ -245,8 +245,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto integrate_planck_from_0_to_x(const // y = x * x * x / 3.0; // 1st order y = (-4 + x) * x + 8 * std::log((2 + x) / 2); // 2nd order // Y_INTERP_MIN is the minimum value returned from interpolate_planck_integral. To ensure y is monotonic with respect to x: + // AMREX_ASSERT(y <= Y_INTERP_MIN); if (y > Y_INTERP_MIN) { y = Y_INTERP_MIN; + } else if (y < 0.) { + y = 0.; } } else if (logx >= LOG_X_MAX) { return 1.0; diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index d0710a6f9..4c6ce87b9 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -110,7 +110,6 @@ template class RadSystem : public HyperbolicSystem::c_hat; static constexpr double radiation_constant_ = RadSystem_Traits::radiation_constant; - static constexpr double Erad_floor_ = RadSystem_Traits::Erad_floor; static constexpr bool compute_v_over_c_terms_ = RadSystem_Traits::compute_v_over_c_terms; static constexpr int nGroups_ = Physics_Traits::nGroups; @@ -122,6 +121,7 @@ template class RadSystem : public HyperbolicSystem::Erad_floor / nGroups_; static constexpr double mean_molecular_mass_ = quokka::EOS_Traits::mean_molecular_mass; static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; @@ -174,6 +174,7 @@ template class RadSystem : public HyperbolicSystem double; AMREX_GPU_HOST_DEVICE static auto ComputePlanckOpacity(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeFluxMeanOpacity(double rho, double Tgas) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto ComputeEnergyMeanOpacity(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputePlanckOpacityTempDerivative(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeEintFromEgas(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Etot) -> double; AMREX_GPU_HOST_DEVICE static auto ComputeEgasFromEint(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Eint) -> double; @@ -186,6 +187,13 @@ template class RadSystem : public HyperbolicSystem const &boundaries, amrex::Real temperature) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray; + + AMREX_GPU_HOST_DEVICE static auto ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) + -> quokka::valarray; + template AMREX_GPU_DEVICE static auto ComputeCellOpticalDepth(const quokka::Array4View &consVar, amrex::GpuArray dx, int i, int j, int k) @@ -218,18 +226,41 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(am radEnergyFractions[g] = y - previous; previous = y; } - // A hack to avoid zeros in radEnergyFractions - // TODO(CCH): use energy_floor - for (int g = 0; g < nGroups_; ++g) { - if (radEnergyFractions[g] < Erad_fraction_floor) { - radEnergyFractions[g] = Erad_fraction_floor; - } - } - radEnergyFractions /= sum(radEnergyFractions); + auto tote = sum(radEnergyFractions); + // AMREX_ALWAYS_ASSERT(tote <= 1.0); + // AMREX_ALWAYS_ASSERT(tote > 0.9999); + radEnergyFractions /= tote; return radEnergyFractions; } } +// define ComputeThermalRadiation, returns the thermal radiation power for each photon group. = a_r * T^4 * radEnergyFractions +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + double power = radiation_constant_ * std::pow(temperature, 4); + auto Erad_g = power * radEnergyFractions; + // set floor + for (int g = 0; g < nGroups_; ++g) { + if (Erad_g[g] < Erad_floor_) { + Erad_g[g] = Erad_floor_; + } + } + return Erad_g; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + // by default, d emission/dT = 4 emission / T + auto erad = ComputeThermalRadiation(temperature, boundaries); + return 4. * erad / temperature; +} + // Linear equation solver for matrix with non-zeros at the first row, first column, and diagonal only. // solve the linear system // [a00 a0i] [x0] = [y0] @@ -919,6 +950,8 @@ void RadSystem::ComputeFluxes(array_t &x1Flux_in, array_t &x1FluxDiff // const quokka::valarray epsilon = {S_corr, S_corr, S_corr, S_corr}; // Jiang et al. (2013) quokka::valarray epsilon = {S_corr * S_corr, S_corr, S_corr, S_corr}; // this code + // fix odd-even instability that appears in the asymptotic diffusion limit + // [for details, see section 3.1: https://ui.adsabs.harvard.edu/abs/2022MNRAS.512.1499R/abstract] if ((i + j + k) % 2 == 1) { // revert to more diffusive flux (has no effect in optically-thin limit) epsilon = {S_corr, 1.0, 1.0, 1.0}; @@ -974,6 +1007,12 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const do return kappaFVec; } +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEnergyMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} + template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double /* rho */, const double /* Tgas */) -> quokka::valarray @@ -1043,7 +1082,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (int g = 0; g < nGroups_; ++g) { Erad0Vec[g] = consPrev(i, j, k, radEnergy_index + numRadVars_ * g); } - AMREX_ASSERT(min(Erad0Vec) >= 0.0); + AMREX_ASSERT(min(Erad0Vec) > 0.0); const double Erad0 = sum(Erad0Vec); // load radiation energy source term @@ -1053,12 +1092,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne Src[g] = dt * (chat * radEnergySource(i, j, k, g)); } - const double a_rad = radiation_constant_; double Egas0 = NAN; double Ekin0 = NAN; double Egas_guess = NAN; double T_gas = NAN; - double fourPiB = NAN; + quokka::valarray fourPiB{}; quokka::valarray EradVec_guess{}; quokka::valarray kappaPVec{}; quokka::valarray kappaEVec{}; @@ -1093,14 +1131,14 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne double deltaEgas = NAN; double Rtot = NAN; double dRtot_dEgas = NAN; - quokka::valarray Rvec; - quokka::valarray dRtot_dErad; - quokka::valarray dRvec_dEgas; - quokka::valarray dFG_dErad; - quokka::valarray dFR_dEgas; - quokka::valarray dFR_i_dErad_i; - quokka::valarray deltaErad; - quokka::valarray F_R; + quokka::valarray Rvec{}; + quokka::valarray dRtot_dErad{}; + quokka::valarray dRvec_dEgas{}; + quokka::valarray dFG_dErad{}; + quokka::valarray dFR_dEgas{}; + quokka::valarray dFR_i_dErad_i{}; + quokka::valarray deltaErad{}; + quokka::valarray F_R{}; const double resid_tol = 1.0e-13; // 1.0e-15; const int maxIter = 400; @@ -1111,30 +1149,29 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(T_gas >= 0.); // compute opacity, emissivity - fourPiB = chat * a_rad * std::pow(T_gas, 4); - auto B = fourPiB / (4.0 * M_PI); + fourPiB = chat * ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); - kappaEVec = kappaFVec; // TODO(CCH): define ComputeEnergyMeanOpacity + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); // compute derivatives w/r/t T_gas - const double dB_dTgas = (4.0 * B) / T_gas; + const auto dfourPiB_dTgas = chat * ComputeThermalRadiationTempDerivative(T_gas, radBoundaries_g_copy); // compute residuals - if constexpr (nGroups_ > 1) { - radEnergyFractions = ComputePlanckEnergyFractions(radBoundaries_g_copy, T_gas); - AMREX_ASSERT(min(radEnergyFractions) > 0.); - } else { - radEnergyFractions[0] = 1.0; - } - - Rvec = dt * rho * (kappaPVec * fourPiB * radEnergyFractions - chat * kappaEVec * EradVec_guess); + const auto absorption = chat * kappaEVec * EradVec_guess; + const auto emission = kappaPVec * fourPiB; + const auto Ediff = emission - absorption; + // auto tot_ediff = sum(abs(Ediff)); + // if ((tot_ediff <= Erad_floor_) || (tot_ediff / (sum(absorption) + sum(emission)) < resid_tol)) { + // break; + // } + Rvec = dt * rho * Ediff; Rtot = sum(Rvec); F_G = Egas_guess - Egas0 + (c / chat) * Rtot; F_R = EradVec_guess - Erad0Vec - (Rvec + Src); - // check if converged - if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * max(abs(F_R)) / Etot0 < resid_tol)) { + // check relative convergence of the residuals + if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * sum(abs(F_R)) / Etot0 < resid_tol)) { break; } @@ -1143,20 +1180,17 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(!std::isnan(sum(dkappaP_dTgas))); // prepare to compute Jacobian elements - const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); + const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T dRtot_dErad = -dt * rho * kappaEVec * chat; - // TODO(CCH): dB_dTgas * radEnergyFractions is a small-dT approximation. Define ComputePlanckEnergyFractionsTempDerivative. Note - // this is accurate in 1-group case. - dRvec_dEgas = dt * rho / c_v * - (4 * M_PI * kappaPVec * dB_dTgas * radEnergyFractions + dkappaP_dTgas * fourPiB * radEnergyFractions - - chat * dkappaE_dTgas * EradVec_guess); - dRtot_dEgas = sum(dRvec_dEgas); - // TODO(CCH): Consider Ben's suggestion of using a global Newton method. + dRvec_dEgas = dt * rho / c_v * (kappaPVec * dfourPiB_dTgas + dkappaP_dTgas * fourPiB - chat * dkappaE_dTgas * EradVec_guess); // Equivalent of eta = (eta > 0.0) ? eta : 0.0, to ensure possitivity of Egas_guess - if (dRtot_dEgas < 0.0) { - dRtot_dEgas = 0.0; + for (int g = 0; g < nGroups_; ++g) { + // AMREX_ASSERT(dRvec_dEgas[g] >= 0.0); + if (dRvec_dEgas[g] < 0.0) { + dRvec_dEgas[g] = 0.0; + } } - // AMREX_ASSERT(dRtot_dEgas >= 0.0); + dRtot_dEgas = sum(dRvec_dEgas); // compute Jacobian elements dFG_dEgas = 1.0 + (c / chat) * dRtot_dEgas; @@ -1173,25 +1207,22 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(!std::isnan(deltaEgas)); AMREX_ASSERT(!deltaErad.hasnan()); - // check relative and absolute convergence of E_r - if ((max(abs(deltaErad / Erad0Vec)) < resid_tol) || (max(abs(deltaErad)) < 0.1 * Erad_floor_)) { - break; - } - EradVec_guess += deltaErad; Egas_guess += deltaEgas; AMREX_ASSERT(min(EradVec_guess) >= 0.); AMREX_ASSERT(Egas_guess > 0.); - // set negative values in EradVec_guess to Erad_floor_ and return the excess to Egas_guess - for (int g = 0; g < nGroups_; ++g) { - if (EradVec_guess[g] < 0.0) { - Egas_guess += (EradVec_guess[g] - Erad_floor_) * (c / chat); - EradVec_guess[g] = Erad_floor_; - } + + // check relative and absolute convergence of E_r + // if ((sum(abs(deltaEgas / Egas_guess)) < 1e-7) || (sum(abs(deltaErad)) <= Erad_floor_)) { + // break; + // } + if (std::abs(deltaEgas / Egas_guess) < 1e-7) { + break; } } // END NEWTON-RAPHSON LOOP AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); + // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); } // endif gamma != 1.0 @@ -1209,24 +1240,17 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::GpuArray Frad_t0{}; amrex::GpuArray Frad_t1{}; amrex::GpuArray dMomentum{0., 0., 0.}; - std::array gasMtm{}; - double realFourPiB = NAN; T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note the T_gas is allowed to be NAN - + quokka::valarray realFourPiB{}; if constexpr (gamma_ != 1.0) { - if constexpr (nGroups_ > 1) { - radEnergyFractions = ComputePlanckEnergyFractions(radBoundaries_g_copy, T_gas); - AMREX_ASSERT(min(radEnergyFractions) > 0.); - } else { - radEnergyFractions[0] = 1.0; - } - kappaPVec = ComputePlanckOpacity(rho, T_gas); - realFourPiB = c * a_rad * std::pow(T_gas, 4); - gasMtm = {x1GasMom0, x2GasMom0, x3GasMom0}; + realFourPiB = c * ComputeThermalRadiation(T_gas, radBoundaries_g_copy); } + kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); + kappaPVec = ComputePlanckOpacity(rho, T_gas); + std::array gasMtm = {x1GasMom0, x2GasMom0, x3GasMom0}; + for (int g = 0; g < nGroups_; ++g) { Frad_t0[0] = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); @@ -1259,7 +1283,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // loop over spatial dimensions for (int n = 0; n < 3; ++n) { - double lastTwoTerms = gasMtm[n] * kappaPVec[g] * realFourPiB * radEnergyFractions[g] * chat / c; + double lastTwoTerms = gasMtm[n] * kappaPVec[g] * realFourPiB[g] * chat / c; // loop over the second rank of the radiation pressure tensor for (int z = 0; z < 3; ++z) { lastTwoTerms += chat * kappaFVec[g] * gasMtm[z] * P[n][z + 1]; @@ -1327,6 +1351,12 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } for (int g = 0; g < nGroups_; ++g) { auto radEnergyNew = EradVec_guess[g] + dErad_work * energyLossFractions[g]; + // AMREX_ASSERT(radEnergyNew > 0.0); + if (radEnergyNew < Erad_floor_) { + // return energy to Egas_guess + Egas_guess -= (Erad_floor_ - radEnergyNew) * (c / chat); + radEnergyNew = Erad_floor_; + } consNew(i, j, k, radEnergy_index + numRadVars_ * g) = radEnergyNew; } } else {