From 66da9804101d281edee00bb0ca010f7d3c832f32 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Tue, 14 May 2024 03:23:39 +1000 Subject: [PATCH] Mean opacity for variable kappa based on piecewise-power-law approximation (#626) ### Description This PR implements piecewise-power-law approximation for multigroup radiation transfer. We calculate Planck-mean, energy-mean, and flux-mean opacity given an arbitrary function kappa(nu, T, rho) based on piecewise-power-law fitting to kappa(nu), E(nu), and F(nu). The user can enable this routine by setting `OpacityModel` to `OpacityModel::piecewisePowerLaw`. See examples in test_radhydro_shock_multigroup.cpp and test_radhydro_pulse_MG_int.hpp. The original approach that the opacities are defined via `ComputePlanckOpacity`, `ComputeFluxMeanOpacity`, etc are retained by setting `OpacityModel` to `OpacityModel::user`, which is the default. We calculate the power-law slope of a radiation quantity via the following relations: $$ \alpha_{i} = {\rm minmod}(s_{i-1/2}, s_{i+1/2}) $$ where $$ s_{i+1/2} = \frac{\ln \frac{X_{i+1}}{\nu_{i+3/2} - \nu_{i+1/2}} - \ln \frac{X_i}{\nu_{i+1/2} - \nu_{i-1/2}}}{\ln \sqrt{\nu_{i+1/2} \nu_{i+3/2}} - \ln \sqrt{\nu_{i-1/2} \nu_{i+1/2}}} $$ I will try to construct a better estimate of the power-law slope in a separate PR. The user can choose to overwrite this by defining your own `RadSystem::ComputeRadQuantityExponents`; see example in test_radhydro_pulse_MG_int.cpp. The power-law fitting of $\kappa(\nu)$ can be specified in two ways. The first method is to specify its power-law slope and lower-bound value for every photon bin, given $T$ and $\rho$. The second method is to specify a precise definition of $\kappa(\nu, T, \rho)$. Then, a power-law fitting to $\kappa(\nu)$ is done on the fly in the following way: $$ \begin{gather} \kappa_{L}^i = \kappa(\nu_{i-1/2}, T, \rho) \\ \alpha^i_{\kappa} = \frac{\ln(\kappa(\nu_{i+1/2}, T, \rho) / \kappa(\nu_{i-1/2},T,\rho))}{\ln(\nu_{i+1/2} / \nu_{i-1/2})} \\ \kappa_{{\rm fit}}^i(\nu) = \kappa_{L}^i \left( \frac{\nu}{\nu_{i-1/2}} \right)^{\alpha_{\kappa}^i} \end{gather} $$ This is definitely not the most accurate way to do this as it is not volume conservative, but it's simple and this is important since it's done for every cell. An alternative is to construct $\kappa_{L,i}$ and $\alpha_{\kappa,i}$ in the beginning if $\kappa$ does not depend on $T$ or $\rho$. ### 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). - [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/RadTube/test_radiation_tube.cpp | 1 + src/RadhydroPulseMGint/CMakeLists.txt | 15 + .../test_radhydro_pulse_MG_int.cpp | 554 ++++++++++++++++++ .../test_radhydro_pulse_MG_int.hpp | 21 + .../test_radhydro_shock_multigroup.cpp | 17 + src/radiation_system.hpp | 335 +++++++++-- 7 files changed, 894 insertions(+), 50 deletions(-) create mode 100644 src/RadhydroPulseMGint/CMakeLists.txt create mode 100644 src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp create mode 100644 src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b08290111..914bdfa6e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -177,6 +177,7 @@ add_subdirectory(RadhydroPulse) add_subdirectory(RadhydroPulseDyn) add_subdirectory(RadhydroPulseGrey) add_subdirectory(RadhydroPulseMG) +add_subdirectory(RadhydroPulseMGint) add_subdirectory(BinaryOrbitCIC) add_subdirectory(Cooling) diff --git a/src/RadTube/test_radiation_tube.cpp b/src/RadTube/test_radiation_tube.cpp index dcc864644..847d63263 100644 --- a/src/RadTube/test_radiation_tube.cpp +++ b/src/RadTube/test_radiation_tube.cpp @@ -65,6 +65,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = C::k_B; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0., 3.3 * T0, inf}; // Kelvin static constexpr int beta_order = 1; + static constexpr OpacityModel opacity_model = OpacityModel::user; }; template <> diff --git a/src/RadhydroPulseMGint/CMakeLists.txt b/src/RadhydroPulseMGint/CMakeLists.txt new file mode 100644 index 000000000..a3efd3bc1 --- /dev/null +++ b/src/RadhydroPulseMGint/CMakeLists.txt @@ -0,0 +1,15 @@ +if (AMReX_SPACEDIM EQUAL 1) + add_executable(test_radhydro_pulse_MG_int test_radhydro_pulse_MG_int.cpp ../fextract.cpp ${QuokkaObjSources}) + + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_pulse_MG_int) + endif(AMReX_GPU_BACKEND MATCHES "CUDA") + + # define variable WORK_PATH = {WORKSPACE}/sims + # set(WORK_PATH "${CMAKE_SOURCE_DIR}/tests/20240408-MG-int-v8.2") + + # mkdir -p WORK_PATH + # file(MAKE_DIRECTORY ${WORK_PATH}) + + add_test(NAME RadhydroPulseMGint COMMAND test_radhydro_pulse_MG_int RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp new file mode 100644 index 000000000..dc2a3b252 --- /dev/null +++ b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp @@ -0,0 +1,554 @@ +/// \file test_radhydro_pulse_MG_int.cpp +/// \brief Defines a test problem for multigroup radiation in the diffusion regime with advection by gas using group-integrated opacity. +/// + +#include "test_radhydro_pulse_MG_int.hpp" +#include "AMReX_BC_TYPES.H" +#include "AMReX_Print.H" +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "physics_info.hpp" +#include "planck_integral.hpp" +#include "radiation_system.hpp" + +struct MGProblem { +}; // dummy type to allow compile-type polymorphism via template specialization +struct ExactProblem { +}; + +// A fixed power law for radiation quantities; for testing purpose only +AMREX_GPU_MANAGED double spec_power = -1.0; // NOLINT +static constexpr bool export_csv = true; + +// constexpr int n_groups_ = 2; +// constexpr amrex::GpuArray rad_boundaries_{1e16, 1e18, 1e20}; +constexpr int n_groups_ = 4; +constexpr amrex::GpuArray rad_boundaries_{1e16, 1e17, 1e18, 1e19, 1e20}; +// constexpr int n_groups_ = 8; +// constexpr amrex::GpuArray rad_boundaries_{1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19, 3.16e19, 1e20}; +// constexpr int n_groups_ = 16; +// constexpr amrex::GpuArray +// rad_boundaries_{1.00000000e+16, 1.77827941e+16, 3.16227766e+16, 5.62341325e+16, 1.00000000e+17, 1.77827941e+17, 3.16227766e+17, 5.62341325e+17, 1.00000000e+18, +// 1.77827941e+18, 3.16227766e+18, 5.62341325e+18, 1.00000000e+19, 1.77827941e+19, 3.16227766e+19, 5.62341325e+19, 1.00000000e+20}; constexpr int n_groups_ = +// 32; constexpr amrex::GpuArray +// rad_boundaries_{1.00000000e+16, 1.33352143e+16, 1.77827941e+16, 2.37137371e+16, 3.16227766e+16, 4.21696503e+16, 5.62341325e+16, 7.49894209e+16, 1.00000000e+17, +// 1.33352143e+17, 1.77827941e+17, 2.37137371e+17, 3.16227766e+17, 4.21696503e+17, 5.62341325e+17, 7.49894209e+17, 1.00000000e+18, 1.33352143e+18, 1.77827941e+18, +// 2.37137371e+18, 3.16227766e+18, 4.21696503e+18, 5.62341325e+18, 7.49894209e+18, 1.00000000e+19, 1.33352143e+19, 1.77827941e+19, 2.37137371e+19, 3.16227766e+19, +// 4.21696503e+19, 5.62341325e+19, 7.49894209e+19, 1.00000000e+20}; constexpr int n_groups_ = 64; constexpr amrex::GpuArray +// rad_boundaries_{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, +// 1.15478198e+19, 1.33352143e+19, 1.53992653e+19, 1.77827941e+19, 2.05352503e+19, 2.37137371e+19, 2.73841963e+19, 3.16227766e+19, 3.65174127e+19, 4.21696503e+19, +// 4.86967525e+19, 5.62341325e+19, 6.49381632e+19, 7.49894209e+19, 8.65964323e+19, 1.00000000e+20}; + +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; + +// static diffusion: (for single group) tau = 2e3, beta = 3e-5, beta tau = 6e-2 +constexpr double kappa0 = 180.; // cm^2 g^-1 +constexpr double v0_adv = 1.0e6; // advecting pulse +constexpr double max_time = 4.8e-6; // max_time = 0.2 * 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); // = 4.799243073 = 1 / 0.2083661912 + +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 = 1; +}; + +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_; + static constexpr int beta_order = 1; + static constexpr OpacityModel opacity_model = OpacityModel::piecewisePowerLaw; +}; +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 int beta_order = 1; + static constexpr OpacityModel opacity_model = OpacityModel::user; +}; + +template <> +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeRadQuantityExponents(ArrayType const & /*quant*/, + amrex::GpuArray const & /*boundaries*/) + -> amrex::GpuArray +{ + amrex::GpuArray exponents{}; + for (int g = 0; g < nGroups_; ++g) { + exponents[g] = spec_power; + } + return exponents; +} + +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_)); +} + +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray const rad_boundaries, const double rho, const double Tgas) + -> amrex::GpuArray, 2> +{ + amrex::GpuArray exponents{}; + amrex::GpuArray kappa_lower{}; + for (int g = 0; g < nGroups_; ++g) { + auto kappa_up = compute_kappa(rad_boundaries[g + 1], Tgas); + auto kappa_down = compute_kappa(rad_boundaries[g], Tgas); + exponents[g] = std::log(kappa_up / kappa_down) / std::log(rad_boundaries[g + 1] / rad_boundaries[g]); + kappa_lower[g] = kappa_down / rho; + AMREX_ASSERT(!std::isnan(exponents[g])); + AMREX_ASSERT(kappa_lower[g] >= 0.); + } + amrex::GpuArray, 2> const exponents_and_values{exponents, kappa_lower}; + return exponents_and_values; +} + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + const double sigma = 3063.96 * std::pow(Tgas / T0, -3.5); + quokka::valarray kappaPVec{}; + kappaPVec.fillin(sigma / rho); + return kappaPVec; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + const double sigma = 101.248 * std::pow(Tgas / T0, -3.5); + quokka::valarray kappaPVec{}; + kappaPVec.fillin(sigma / rho); + return kappaPVec; +} + +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 + 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]); + + // 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 Erad = a_rad * std::pow(Trad, 4); + const double rho = compute_exact_rho(x - x0); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double v0 = v0_adv; + + // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; + state_cc(i, j, k, RadSystem::radEnergy_index) = Erad; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = 4. / 3. * v0 * Erad; + 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) = 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 based on the radhydro_pulse test and is a test of interpolation for variable opacity for multigroup radiation. + + // 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 + + amrex::ParmParse const pp("rad"); + pp.query("spec_power", spec_power); + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + for (int i = 0; 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: advecting pulse with multigroup integration + + // 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); + int nx = static_cast(position.size()); + amrex::GpuArray prob_lo = sim.geom[0].ProbLoArray(); + amrex::GpuArray prob_hi = sim.geom[0].ProbHiArray(); + // compute the pixel size + double dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx); + double move = v0_adv * sim.tNew_[0]; + int n_p = static_cast(move / dx); + int half = static_cast(nx / 2.0); + double drift = move - static_cast(n_p) * dx; + int shift = n_p - static_cast((n_p + half) / nx) * nx; + + 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) { + 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); + } + } + amrex::Real const x = position[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]; + xs.at(i) = x - drift; + rhogas.at(index_) = rho_t; + Trad.at(index_) = Trad_t; + Tgas.at(index_) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); + Vgas.at(index_) = 1e-5 * (v_t - v0_adv); + } + // END OF PROBLEM 1 + + // Problem 2: exact opacity + + // 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(); + + auto [position0, values0] = fextract(sim2.state_new_cc_[0], sim2.Geom(0), 0, 0.0); + + // evolve + sim2.evolve(); + + // read output variables + auto [position2, values2] = fextract(sim2.state_new_cc_[0], sim2.Geom(0), 0, 0.0); + nx = static_cast(position2.size()); + prob_lo = sim2.geom[0].ProbLoArray(); + prob_hi = sim2.geom[0].ProbHiArray(); + // compute the pixel size + dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx); + move = v0_adv * sim2.tNew_[0]; + n_p = static_cast(move / dx); + half = static_cast(nx / 2.0); + drift = move - static_cast(n_p) * dx; + shift = n_p - static_cast((n_p + half) / nx) * nx; + + std::vector xs2(nx); + std::vector Trad2(nx); + std::vector xs0(nx); + std::vector Trad0(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]; + const auto Erad_0 = values0.at(RadSystem::radEnergy_index)[i]; + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + const auto Trad_0 = std::pow(Erad_0 / 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); + const auto x0 = position0[i]; + xs0.at(i) = x0; + Trad0.at(i) = Trad_0; + } + // 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] - Tgas2[i]); + sol_norm += std::abs(Tgas2[i]); + } + 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 (MG integrated)"; + Trad_args["linestyle"] = "-"; + Trad_args["color"] = "C0"; + Tgas_args["label"] = "Tgas (MG integrated)"; + Tgas_args["linestyle"] = "--"; + Tgas_args["color"] = "C1"; + matplotlibcpp::plot(xs, Trad, Trad_args); + matplotlibcpp::plot(xs, Tgas, Tgas_args); + Trad_args["label"] = "Trad (grey)"; + Trad_args["linestyle"] = "-"; + Trad_args["color"] = "k"; + Tgas_args["label"] = "Tgas (grey)"; + Tgas_args["linestyle"] = "--"; + Tgas_args["color"] = "grey"; + matplotlibcpp::plot(xs2, Trad2, Trad_args); + matplotlibcpp::plot(xs2, Tgas2, Tgas_args); + // matplotlibcpp::ylim(0.98e7, 1.35e7); + matplotlibcpp::ylim(0.98e7, 2.02e7); + matplotlibcpp::grid(true); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("temperature (K)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("nGroups = {}, time t = {:.4g}", n_groups_, sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_int_temperature.pdf"); + + // plot gas density profile + matplotlibcpp::clf(); + std::map rho_args; + rho_args["label"] = "MG integrated"; + rho_args["linestyle"] = "-"; + rho_args["color"] = "C0"; + matplotlibcpp::plot(xs, rhogas, rho_args); + rho_args["label"] = "grey)"; + rho_args["linestyle"] = "--"; + rho_args["color"] = "k"; + matplotlibcpp::plot(xs2, rhogas2, rho_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("gas density (g cm^-3)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("nGroups = {}, time t = {:.4g}", n_groups_, sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_int_density.pdf"); + + // plot gas velocity profile + matplotlibcpp::clf(); + std::map vgas_args; + vgas_args["label"] = "MG interpolated"; + vgas_args["linestyle"] = "-"; + vgas_args["color"] = "C0"; + matplotlibcpp::plot(xs, Vgas, vgas_args); + vgas_args["label"] = "grey"; + vgas_args["linestyle"] = "--"; + vgas_args["color"] = "k"; + matplotlibcpp::plot(xs2, Vgas2, vgas_args); + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("gas velocity (km s^-1)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("nGroups = {}, time t = {:.4g}", n_groups_, sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radhydro_pulse_MG_int_velocity.pdf"); +#endif + + // Save xs, Trad, Tgas, rhogas, Vgas, xs_mg, Trad_mg, Tgas_mg, rhogas_mg, Vgas_mg, xs2, Trad2, Tgas2, rhogas2, Vgas2 + if (export_csv) { + std::ofstream file; + file.open("radhydro_pulse_MG_int_temperature.csv"); + file << "xs,Trad,Tgas,xs2,Trad2,Tgas2\n"; + for (size_t i = 0; i < xs.size(); ++i) { + file << std::scientific << std::setprecision(12) << xs[i] << "," << Trad[i] << "," << Tgas[i] << "," << xs2[i] << "," << Trad2[i] << "," + << Tgas2[i] << "\n"; + } + file.close(); + + // Save xs, rhogas, xs_mg, rhogas_mg, xs2, rhogas2 + file.open("radhydro_pulse_MG_int_density.csv"); + file << "xs,rhogas,xs2,rhogas2\n"; + for (size_t i = 0; i < xs.size(); ++i) { + file << std::scientific << std::setprecision(12) << xs[i] << "," << rhogas[i] << "," << xs2[i] << "," << rhogas2[i] << "\n"; + } + file.close(); + + // Save xs, Vgas, xs_mg, Vgas_mg, xs2, Vgas2 + file.open("radhydro_pulse_MG_int_velocity.csv"); + file << "xs,Vgas,xs2,Vgas2\n"; + for (size_t i = 0; i < xs.size(); ++i) { + file << std::scientific << std::setprecision(12) << xs[i] << "," << Vgas[i] << "," << xs2[i] << "," << Vgas2[i] << "\n"; + } + file.close(); + } + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + return status; +} diff --git a/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.hpp b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.hpp new file mode 100644 index 000000000..a077f921a --- /dev/null +++ b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.hpp @@ -0,0 +1,21 @@ +#ifndef TEST_RADHYDRO_PULSE_MG_INT_HPP_ // NOLINT +#define TEST_RADHYDRO_PULSE_MG_INT_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_INT_HPP_ diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index 446be45f8..0fe0cacbf 100644 --- a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -9,6 +9,7 @@ #include "ArrayUtil.hpp" #include "fextract.hpp" +#include "radiation_system.hpp" #include "test_radhydro_shock_multigroup.hpp" struct ShockProblem { @@ -62,6 +63,7 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.00000000e+15, 1.00000000e+16, 1.00000000e+17, 1.00000000e+18, 1.00000000e+19, 1.00000000e+20}; static constexpr int beta_order = 1; + static constexpr OpacityModel opacity_model = OpacityModel::piecewisePowerLaw; }; template <> struct quokka::EOS_Traits { @@ -70,6 +72,21 @@ template <> struct quokka::EOS_Traits { static constexpr double gamma = gamma_gas; }; +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, const double /*Tgas*/) + -> amrex::GpuArray, 2> +{ + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_; ++i) { + exponents_and_values[0][i] = 0.0; + } + for (int i = 0; i < nGroups_; ++i) { + exponents_and_values[1][i] = kappa / rho; + } + return exponents_and_values; +} + template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> quokka::valarray diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index ea7ba04d8..b3fb6e27b 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -19,6 +19,7 @@ #include "AMReX_Array.H" #include "AMReX_BLassert.H" #include "AMReX_GpuQualifiers.H" +#include "AMReX_IParser_Y.H" #include "AMReX_REAL.H" // internal headers @@ -50,6 +51,13 @@ static constexpr double inf = std::numeric_limits::max(); // Optional: include a wavespeed correction term in the radiation flux to suppress instability static const bool use_wavespeed_correction = false; +// enum for opacity_model +enum class OpacityModel { + user = 0, // user-defined opacity for each group, given as a function of density and temperature. + piecewisePowerLaw // piecewise power-law opacity model with piecewise power-law fitting to a user-defined opacity function and on-the-fly piecewise + // power-law fitting to radiation energy density and flux. +}; + // this struct is specialized by the user application code // template struct RadSystem_Traits { @@ -60,6 +68,7 @@ template struct RadSystem_Traits { static constexpr double energy_unit = C::ev2erg; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries = {0., inf}; static constexpr double beta_order = 1; + static constexpr OpacityModel opacity_model = OpacityModel::user; }; // A struct to hold the results of the ComputeRadPressure function. @@ -68,6 +77,19 @@ struct RadPressureResult { double S; // maximum wavespeed for the radiation system }; +[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto minmod_func(double a, double b) -> double +{ + return 0.5 * (sgn(a) + sgn(b)) * std::min(std::abs(a), std::abs(b)); +} + +// Use SFINAE (Substitution Failure Is Not An Error) to check if opacity_model is defined in RadSystem_Traits +template struct RadSystem_Has_Opacity_Model : std::false_type { +}; + +template +struct RadSystem_Has_Opacity_Model::opacity_model)>> : std::true_type { +}; + /// Class for the radiation moment equations /// template class RadSystem : public HyperbolicSystem @@ -122,6 +144,14 @@ template class RadSystem : public HyperbolicSystem::Erad_floor / nGroups_; + static constexpr OpacityModel opacity_model_ = []() constexpr { + if constexpr (RadSystem_Has_Opacity_Model::value) { + return RadSystem_Traits::opacity_model; + } else { + return OpacityModel::user; + } + }(); + static constexpr double mean_molecular_mass_ = quokka::EOS_Traits::mean_molecular_mass; static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; static constexpr double gamma_ = quokka::EOS_Traits::gamma; @@ -164,17 +194,26 @@ template class RadSystem : public HyperbolicSystem 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 DefineOpacityExponentsAndLowerValues(amrex::GpuArray rad_boundaries, double rho, double Tgas) + -> amrex::GpuArray, 2>; + AMREX_GPU_HOST_DEVICE static auto ComputeGroupMeanOpacity(amrex::GpuArray, 2> kappa_expo_and_lower_value, + amrex::GpuArray radBoundaryRatios, + amrex::GpuArray alpha_quant) -> 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; + template + AMREX_GPU_HOST_DEVICE static auto ComputeRadQuantityExponents(ArrayType const &quant, amrex::GpuArray const &boundaries) + -> amrex::GpuArray; + AMREX_GPU_HOST_DEVICE static void SolveLinearEqs(double a00, const quokka::valarray &a0i, const quokka::valarray &ai0, const quokka::valarray &aii, const double &y0, const quokka::valarray &yi, double &x0, quokka::valarray &xi); AMREX_GPU_HOST_DEVICE static auto Solve3x3matrix(double C00, double C01, double C02, double C10, double C11, double C12, double C20, double C21, - double C22, double Y0, double Y1, double Y2) -> quokka::valarray; + double C22, double Y0, double Y1, double Y2) -> std::tuple; AMREX_GPU_HOST_DEVICE static auto ComputePlanckEnergyFractions(amrex::GpuArray const &boundaries, amrex::Real temperature) -> quokka::valarray; @@ -274,12 +313,11 @@ AMREX_GPU_HOST_DEVICE void RadSystem::SolveLinearEqs(const double a00 template AMREX_GPU_HOST_DEVICE auto RadSystem::Solve3x3matrix(const double C00, const double C01, const double C02, const double C10, const double C11, const double C12, const double C20, const double C21, const double C22, const double Y0, - const double Y1, const double Y2) -> quokka::valarray + const double Y1, const double Y2) -> std::tuple { // Solve the 3x3 matrix equation: C * X = Y under the assumption that only the diagonal terms // are guaranteed to be non-zero and are thus allowed to be divided by. - quokka::valarray ret{}; auto E11 = C11 - C01 * C10 / C00; auto E12 = C12 - C02 * C10 / C00; auto E21 = C21 - C01 * C20 / C00; @@ -289,8 +327,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::Solve3x3matrix(const double C00 auto X2 = (Z2 - Z1 * E21 / E11) / (E22 - E12 * E21 / E11); auto X1 = (Z1 - E12 * X2) / E11; auto X0 = (Y0 - C01 * X1 - C02 * X2) / C00; - ret = {X0, X1, X2}; - return ret; + + return std::make_tuple(X0, X1, X2); } template @@ -861,6 +899,89 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEnergyMeanOpacity(const return ComputePlanckOpacity(rho, Tgas); } +template +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, + const double /*rho*/, const double /*Tgas*/) + -> amrex::GpuArray, 2> +{ + amrex::GpuArray, 2> exponents_and_values{}; + return exponents_and_values; +} + +template +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeRadQuantityExponents(ArrayType const &quant, amrex::GpuArray const &boundaries) + -> amrex::GpuArray +{ + // Compute the exponents for the radiation energy density, radiation flux, radiation pressure, or Planck function. + + // Note: Could save some memory by using bin_center_previous and bin_center_current + amrex::GpuArray bin_center{}; + amrex::GpuArray quant_mean{}; + amrex::GpuArray logslopes{}; + amrex::GpuArray exponents{}; + for (int g = 0; g < nGroups_; ++g) { + bin_center[g] = std::sqrt(boundaries[g] * boundaries[g + 1]); + quant_mean[g] = quant[g] / (boundaries[g + 1] - boundaries[g]); + if (g > 0) { + AMREX_ASSERT(bin_center[g] > bin_center[g - 1]); + if (quant_mean[g] == 0.0 && quant_mean[g - 1] == 0.0) { + logslopes[g - 1] = 0.0; + } else if (quant_mean[g - 1] * quant_mean[g] <= 0.0) { + if (quant_mean[g] > quant_mean[g - 1]) { + logslopes[g - 1] = inf; + } else { + logslopes[g - 1] = -inf; + } + } else { + logslopes[g - 1] = std::log(std::abs(quant_mean[g] / quant_mean[g - 1])) / std::log(bin_center[g] / bin_center[g - 1]); + } + AMREX_ASSERT(!std::isnan(logslopes[g - 1])); + } + } + for (int g = 0; g < nGroups_; ++g) { + if (g == 0 || g == nGroups_ - 1) { + exponents[g] = 0.0; + } else { + exponents[g] = minmod_func(logslopes[g - 1], logslopes[g]); + } + AMREX_ASSERT(!std::isnan(exponents[g])); + AMREX_ASSERT(std::abs(exponents[g]) < 100); + } + + return exponents; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeGroupMeanOpacity(amrex::GpuArray, 2> kappa_expo_and_lower_value, + amrex::GpuArray radBoundaryRatios, + amrex::GpuArray alpha_quant) -> quokka::valarray +{ + amrex::GpuArray const &alpha_kappa = kappa_expo_and_lower_value[0]; + amrex::GpuArray const &kappa_lower = kappa_expo_and_lower_value[1]; + + quokka::valarray kappa{}; + for (int g = 0; g < nGroups_; ++g) { + double alpha = alpha_quant[g] + 1.0; + double part1 = 0.0; + if (std::abs(alpha) < 1e-8) { + part1 = std::log(radBoundaryRatios[g]); + } else { + part1 = (std::pow(radBoundaryRatios[g], alpha) - 1.0) / alpha; + } + alpha = alpha_quant[g] + alpha_kappa[g] + 1.0; + double part2 = 0.0; + if (std::abs(alpha) < 1e-8) { + part2 = std::log(radBoundaryRatios[g]); + } else { + part2 = (std::pow(radBoundaryRatios[g], alpha) - 1.0) / alpha; + } + kappa[g] = kappa_lower[g] / part1 * part2; + AMREX_ASSERT(!std::isnan(kappa[g])); + } + return kappa; +} + template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double /* rho */, const double /* Tgas */) -> quokka::valarray @@ -902,9 +1023,14 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne dt = (1.0 - IMEX_a32) * dt_radiation; } - amrex::GpuArray radBoundaries_g{}; + amrex::GpuArray radBoundaries_g = radBoundaries_; + amrex::GpuArray radBoundaryRatios{}; if constexpr (nGroups_ > 1) { - radBoundaries_g = RadSystem_Traits::radBoundaries; + if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + for (int g = 0; g < nGroups_; ++g) { + radBoundaryRatios[g] = radBoundaries_g[g + 1] / radBoundaries_g[g]; + } + } } // Add source terms @@ -954,12 +1080,19 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray kappaPVec{}; quokka::valarray kappaEVec{}; quokka::valarray kappaFVec{}; + amrex::GpuArray, 2> kappa_expo_and_lower_value{}; + amrex::GpuArray alpha_B{}; + amrex::GpuArray alpha_E{}; + amrex::GpuArray alpha_F{}; quokka::valarray kappaPoverE{}; quokka::valarray tau0{}; // optical depth across c * dt at old state quokka::valarray tau{}; // optical depth across c * dt at new state quokka::valarray D{}; // D = S / tau0 quokka::valarray work{}; quokka::valarray work_prev{}; + amrex::GpuArray, 3> frad{}; + amrex::GpuArray dMomentum{}; + amrex::GpuArray, 3> Frad_t1{}; work.fillin(0.0); work_prev.fillin(0.0); @@ -973,13 +1106,14 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // make a copy of radBoundaries_g amrex::GpuArray radBoundaries_g_copy{}; + amrex::GpuArray radBoundaryRatios_copy{}; for (int g = 0; g < nGroups_ + 1; ++g) { radBoundaries_g_copy[g] = radBoundaries_g[g]; + if (g < nGroups_) { + radBoundaryRatios_copy[g] = radBoundaryRatios[g]; + } } - amrex::GpuArray dMomentum{}; - amrex::GpuArray, nGroups_> Frad_t1{}; - amrex::Real gas_update_factor = 1.0; if (stage == 1) { gas_update_factor = IMEX_a32; @@ -1039,11 +1173,21 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); AMREX_ASSERT(T_gas >= 0.); fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); + + if constexpr (opacity_model_ == OpacityModel::user) { + kappaPVec = ComputePlanckOpacity(rho, T_gas); + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); + kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + } AMREX_ASSERT(!kappaPVec.hasnan()); AMREX_ASSERT(!kappaEVec.hasnan()); + AMREX_ASSERT(!kappaFVec.hasnan()); for (int g = 0; g < nGroups_; ++g) { if (kappaEVec[g] > 0.0) { @@ -1057,13 +1201,36 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // compute the work term at the old state // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); if (ite == 0) { - for (int g = 0; g < nGroups_; ++g) { - // work[g] = dt * chat * rho * kappaPVec[g] * (Erad0Vec[g] - fourPiBoverC[g]); - const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); - const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); - const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); - work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * chat / (c * c) * - lorentz_factor_v * kappaFVec[g] * dt; + if constexpr (opacity_model_ == OpacityModel::user) { + for (int g = 0; g < nGroups_; ++g) { + // work[g] = dt * chat * rho * kappaPVec[g] * (Erad0Vec[g] - fourPiBoverC[g]); + const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + // work = v * F * chi + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * + (2.0 * kappaEVec[g] - kappaFVec[g]); + work[g] *= chat / (c * c) * lorentz_factor_v * dt; + } + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + for (int g = 0; g < nGroups_; ++g) { + frad[0][g] = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + frad[1][g] = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + frad[2][g] = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + work[g] = 0.0; + } + for (int n = 0; n < 3; ++n) { + alpha_F = ComputeRadQuantityExponents(frad[n], radBoundaries_g_copy); + kappaFVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); + for (int g = 0; g < nGroups_; ++g) { + work[g] += + (kappa_expo_and_lower_value[0][g] + 1.0) * gasMtm0[n] * kappaFVec[g] * frad[n][g]; + } + } + for (int g = 0; g < nGroups_; ++g) { + work[g] *= chat / (c * c) * dt; + } } } } @@ -1089,7 +1256,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray deltaD{}; quokka::valarray F_D{}; - const double resid_tol = 1.0e-13; // 1.0e-15; + const double resid_tol = 1.0e-11; // 1.0e-15; const int maxIter = 400; int n = 0; for (; n < maxIter; ++n) { @@ -1098,9 +1265,20 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(T_gas >= 0.); // compute opacity, emissivity fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); + + if constexpr (opacity_model_ == OpacityModel::user) { + kappaPVec = ComputePlanckOpacity(rho, T_gas); + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + } AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); + for (int g = 0; g < nGroups_; ++g) { if (kappaEVec[g] > 0.0) { kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; @@ -1119,11 +1297,19 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne EradVec_guess[g] = kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work[g]) / tau[g]); } } - F_G = Egas_guess - Egas0 + (c / chat) * sum(Rvec); + // F_G = Egas_guess - Egas0 + (c / chat) * sum(Rvec); + F_G = Egas_guess - Egas0; F_D = EradVec_guess - Erad0Vec - (Rvec + Src); + double F_D_abs_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] > 0.0) { + F_D_abs_sum += std::abs(F_D[g]); + F_G += (c / chat) * Rvec[g]; + } + } // check relative convergence of the residuals - if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * sum(abs(F_D)) / Etot0 < resid_tol)) { + if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * F_D_abs_sum / Etot0 < resid_tol)) { break; } @@ -1185,10 +1371,30 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); if constexpr (gamma_ != 1.0) { fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); } - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note that kappaFVec is used no matter what the value of gamma is + + if constexpr (opacity_model_ == OpacityModel::user) { + if constexpr (gamma_ != 1.0) { + kappaPVec = ComputePlanckOpacity(rho, T_gas); + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); + } + kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note that kappaFVec is used no matter what the value of gamma is + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + if constexpr (gamma_ != 1.0) { + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); + } + // Note that alpha_F has not been changed in the Newton iteration + kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); + } + AMREX_ASSERT(!kappaFVec.hasnan()); dMomentum = {0., 0., 0.}; @@ -1214,11 +1420,16 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (int n = 0; n < 3; ++n) { // compute thermal radiation term - double v_term = kappaPVec[g] * fourPiBoverC[g] * lorentz_factor_v; + double v_term = NAN; - // compute (kappa_F - kappa_E) term - if (kappaFVec[g] != kappaEVec[g]) { - v_term += (kappaFVec[g] - kappaEVec[g]) * erad * std::pow(lorentz_factor_v, 3); + if constexpr (opacity_model_ == OpacityModel::user) { + v_term = kappaPVec[g] * fourPiBoverC[g] * lorentz_factor_v; + // compute (kappa_F - kappa_E) term + if (kappaFVec[g] != kappaEVec[g]) { + v_term += (kappaFVec[g] - kappaEVec[g]) * erad * std::pow(lorentz_factor_v, 3); + } + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + v_term = kappaPVec[g] * fourPiBoverC[g] * (2.0 - kappa_expo_and_lower_value[0][g] - alpha_B[g]) / 3.0; } v_term *= chat * dt * gasMtm0[n]; @@ -1228,7 +1439,12 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (int z = 0; z < 3; ++z) { pressure_term += gasMtm0[z] * Tedd[n][z] * erad; } - pressure_term *= chat * dt * kappaFVec[g] * lorentz_factor_v; + if constexpr (opacity_model_ == OpacityModel::user) { + pressure_term *= chat * dt * kappaFVec[g] * lorentz_factor_v; + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + pressure_term *= chat * dt * kappaEVec[g] * (kappa_expo_and_lower_value[0][g] + 1.0); + } + v_term += pressure_term; v_terms[n] = v_term; } @@ -1236,19 +1452,19 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne if constexpr (beta_order_ == 1) { for (int n = 0; n < 3; ++n) { // Compute flux update - Frad_t1[g][n] = (Frad_t0[n] + v_terms[n]) / (1.0 + F_coeff); + Frad_t1[n][g] = (Frad_t0[n] + v_terms[n]) / (1.0 + F_coeff); // Compute conservative gas momentum update - dMomentum[n] += -(Frad_t1[g][n] - Frad_t0[n]) / (c * chat); + dMomentum[n] += -(Frad_t1[n][g] - Frad_t0[n]) / (c * chat); } } else { if (kappaFVec[g] == kappaEVec[g]) { for (int n = 0; n < 3; ++n) { // Compute flux update - Frad_t1[g][n] = (Frad_t0[n] + v_terms[n]) / (1.0 + F_coeff); + Frad_t1[n][g] = (Frad_t0[n] + v_terms[n]) / (1.0 + F_coeff); // Compute conservative gas momentum update - dMomentum[n] += -(Frad_t1[g][n] - Frad_t0[n]) / (c * chat); + dMomentum[n] += -(Frad_t1[n][g] - Frad_t0[n]) / (c * chat); } } else { const double K0 = @@ -1276,24 +1492,27 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const double B1 = v_terms[1] + Frad_t0[1]; const double B2 = v_terms[2] + Frad_t0[2]; - Frad_t1[g] = Solve3x3matrix(A00, A01, A02, A10, A11, A12, A20, A21, A22, B0, B1, B2); + auto [sol0, sol1, sol2] = Solve3x3matrix(A00, A01, A02, A10, A11, A12, A20, A21, A22, B0, B1, B2); + Frad_t1[0][g] = sol0; + Frad_t1[1][g] = sol1; + Frad_t1[2][g] = sol2; for (int n = 0; n < 3; ++n) { - dMomentum[n] += -(Frad_t1[g][n] - Frad_t0[n]) / (c * chat); + dMomentum[n] += -(Frad_t1[n][g] - Frad_t0[n]) / (c * chat); } } } } else { for (int n = 0; n < 3; ++n) { - Frad_t1[g][n] = Frad_t0[n] / (1.0 + rho * kappaFVec[g] * chat * dt); + Frad_t1[n][g] = Frad_t0[n] / (1.0 + rho * kappaFVec[g] * chat * dt); // Compute conservative gas momentum update - dMomentum[n] += -(Frad_t1[g][n] - Frad_t0[n]) / (c * chat); + dMomentum[n] += -(Frad_t1[n][g] - Frad_t0[n]) / (c * chat); } } } // end loop over radiation groups for flux update - amrex::Real x1GasMom1 = consPrev(i, j, k, x1GasMomentum_index) + dMomentum[0]; - amrex::Real x2GasMom1 = consPrev(i, j, k, x2GasMomentum_index) + dMomentum[1]; - amrex::Real x3GasMom1 = consPrev(i, j, k, x3GasMomentum_index) + dMomentum[2]; + amrex::Real const x1GasMom1 = consPrev(i, j, k, x1GasMomentum_index) + dMomentum[0]; + amrex::Real const x2GasMom1 = consPrev(i, j, k, x2GasMomentum_index) + dMomentum[1]; + amrex::Real const x3GasMom1 = consPrev(i, j, k, x3GasMomentum_index) + dMomentum[2]; // 3. Deal with the work term. if constexpr ((gamma_ != 1.0) && (beta_order_ != 0)) { @@ -1320,7 +1539,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // compute energyLossFractions for (int g = 0; g < nGroups_; ++g) { energyLossFractions[g] = - kappaFVec[g] * (x1GasMom1 * Frad_t1[g][0] + x2GasMom1 * Frad_t1[g][1] + x3GasMom1 * Frad_t1[g][2]); + kappaFVec[g] * (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]); } auto energyLossFractionsTot = sum(energyLossFractions); if (energyLossFractionsTot != 0.0) { @@ -1352,8 +1571,24 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // copy work to work_prev work_prev[g] = work[g]; // compute new work term from the updated radiation flux and velocity - work[g] = (x1GasMom1 * Frad_t1[g][0] + x2GasMom1 * Frad_t1[g][1] + x3GasMom1 * Frad_t1[g][2]) * chat / (c * c) * - lorentz_factor_v * kappaFVec[g] * dt; + if constexpr (opacity_model_ == OpacityModel::user) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * chat / (c * c) * + lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; + } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + for (int n = 0; n < 3; ++n) { + work[n] = 0.0; + } + for (int n = 0; n < 3; ++n) { + alpha_F = ComputeRadQuantityExponents(Frad_t1[n], radBoundaries_g_copy); + kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); + for (int g = 0; g < nGroups_; ++g) { + work[g] += (kappa_expo_and_lower_value[0][g] + 1.0) * gasMtm0[n] * kappaFVec[g] * Frad_t1[n][g]; + } + } + for (int g = 0; g < nGroups_; ++g) { + work[g] *= chat * dt / (c * c); + } + } } // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop @@ -1387,9 +1622,9 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne if constexpr (gamma_ != 1.0) { consNew(i, j, k, radEnergy_index + numRadVars_ * g) = EradVec_guess[g]; } - consNew(i, j, k, x1RadFlux_index + numRadVars_ * g) = Frad_t1[g][0]; - consNew(i, j, k, x2RadFlux_index + numRadVars_ * g) = Frad_t1[g][1]; - consNew(i, j, k, x3RadFlux_index + numRadVars_ * g) = Frad_t1[g][2]; + consNew(i, j, k, x1RadFlux_index + numRadVars_ * g) = Frad_t1[0][g]; + consNew(i, j, k, x2RadFlux_index + numRadVars_ * g) = Frad_t1[1][g]; + consNew(i, j, k, x3RadFlux_index + numRadVars_ * g) = Frad_t1[2][g]; } }); }