diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index a16bca36a..601915334 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -1617,10 +1617,10 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real // new radiation state is stored in state_new_cc_ // new hydro state is stored in state_new_cc_ (always the case during radiation update) - amrex::Gpu::Buffer iteration_failure_counter( - {0, 0, 0}); // failure counter for: matter-radiation coupling, dust temperature, outer iteration - amrex::Gpu::Buffer iteration_counter( - {0, 0, 0}); // iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations + // failure counter for: matter-radiation coupling, dust temperature, outer iteration + amrex::Gpu::Buffer iteration_failure_counter({0, 0, 0}); + // iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations, decoupled gas-dust update + amrex::Gpu::Buffer iteration_counter({0, 0, 0, 0}); int *p_iteration_failure_counter = iteration_failure_counter.data(); int *p_iteration_counter = iteration_counter.data(); @@ -1659,13 +1659,15 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real if (print_rad_counter_) { auto h_iteration_counter = iteration_counter.copyToHost(); - long global_solver_count = h_iteration_counter[0]; // number of Newton-Raphson solvings - long global_iteration_sum = h_iteration_counter[1]; // sum of Newton-Raphson iterations - int global_iteration_max = h_iteration_counter[2]; // max number of Newton-Raphson iterations + long global_solver_count = h_iteration_counter[0]; // number of Newton-Raphson solvings + long global_iteration_sum = h_iteration_counter[1]; // sum of Newton-Raphson iterations + int global_iteration_max = h_iteration_counter[2]; // max number of Newton-Raphson iterations + long global_decoupled_iteration_sum = h_iteration_counter[3]; // sum of decoupled gas-dust Newton-Raphson iterations amrex::ParallelDescriptor::ReduceLongSum(global_solver_count); amrex::ParallelDescriptor::ReduceLongSum(global_iteration_sum); amrex::ParallelDescriptor::ReduceIntMax(global_iteration_max); + amrex::ParallelDescriptor::ReduceLongSum(global_decoupled_iteration_sum); if (amrex::ParallelDescriptor::IOProcessor()) { const auto n_cells = CountCells(lev); @@ -1675,9 +1677,15 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real static_cast(global_iteration_sum) / static_cast(global_solver_count); const double global_solving_mean = static_cast(global_solver_count) / static_cast(n_cells) / 2.0; // 2 stages - amrex::Print() << "average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean + const double global_decoupled_iteration_mean = + static_cast(global_decoupled_iteration_sum) / static_cast(global_solver_count); + amrex::Print() << "The average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean << ", (mean, max) number of Newton-Raphson iterations are " << global_iteration_mean << ", " - << global_iteration_max << "\n"; + << global_iteration_max << ".\n"; + if constexpr (RadSystem_Traits::enable_dust_gas_thermal_coupling_model) { + amrex::Print() << "The fraction of gas-dust interactions that are decoupled is " + << global_decoupled_iteration_mean << "\n"; + } } } } diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index 2c6e32528..d2d4195b0 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -36,6 +36,7 @@ add_subdirectory(RadTube) add_subdirectory(RadDust) add_subdirectory(RadDustMG) add_subdirectory(RadMarshakDust) +add_subdirectory(RadMarshakDustPE) add_subdirectory(RadhydroShell) add_subdirectory(RadhydroShock) diff --git a/src/problems/RadDustMG/test_rad_dust_MG.hpp b/src/problems/RadDustMG/test_rad_dust_MG.hpp index 6187b8b7f..bcd899616 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.hpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.hpp @@ -14,7 +14,7 @@ // internal headers #include "math/interpolate.hpp" -#include "radiation/radiation_system.hpp" +#include "radiation/radiation_dust_system.hpp" // function definitions diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp index bd3af42bf..8f39445bb 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -13,7 +13,7 @@ #include "util/fextract.hpp" #include "util/valarray.hpp" -struct StreamingProblem { +struct MarshakProblem { }; AMREX_GPU_MANAGED double kappa1 = NAN; // dust opacity at IR @@ -42,13 +42,13 @@ constexpr int n_group_ = 2; static constexpr amrex::GpuArray radBoundaries_{1e-10, 100, 1e4}; static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity; -template <> struct quokka::EOS_Traits { +template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; -template <> struct Physics_Traits { +template <> struct Physics_Traits { // cell-centred static constexpr bool is_hydro_enabled = false; static constexpr int numMassScalars = 0; // number of mass scalars @@ -59,7 +59,7 @@ template <> struct Physics_Traits { static constexpr int nGroups = n_group_; // number of radiation groups }; -template <> struct RadSystem_Traits { +template <> struct RadSystem_Traits { static constexpr double c_light = c; static constexpr double c_hat = chat; static constexpr double radiation_constant = a_rad; @@ -71,20 +71,25 @@ template <> struct RadSystem_Traits { static constexpr OpacityModel opacity_model = opacity_model_; }; -template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real +template <> struct ISM_Traits { + static constexpr double gas_dust_coupling_threshold = 1.0e-5; + static constexpr bool enable_photoelectric_heating = false; +}; + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return kappa1; } -template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return kappa1; } template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto -RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double /*rho*/, - const double /*Tgas*/) -> amrex::GpuArray, 2> +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_ + 1; ++i) { @@ -98,36 +103,36 @@ RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArra return exponents_and_values; } -template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) { const amrex::Box &indexRange = grid_elem.indexRange_; const amrex::Array4 &state_cc = grid_elem.array_; const auto Egas0 = initial_T * CV; - const auto Erads = RadSystem::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_); + const auto Erads = RadSystem::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_); // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g]; - state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; - 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; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + 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) = Egas0; - state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; - state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas0; - state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; - state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; - state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas0; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; }); } template <> AMREX_GPU_DEVICE AMREX_FORCE_INLINE void -AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, int /*dcomp*/, - int /*numcomp*/, amrex::GeometryData const &geom, const amrex::Real /*time*/, - const amrex::BCRec * /*bcr*/, int /*bcomp*/, int /*orig_comp*/) +AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, int /*dcomp*/, int /*numcomp*/, + amrex::GeometryData const &geom, const amrex::Real /*time*/, const amrex::BCRec * /*bcr*/, + int /*bcomp*/, int /*orig_comp*/) { #if (AMREX_SPACEDIM == 1) auto i = iv.toArray()[0]; @@ -145,7 +150,7 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVec amrex::Box const &box = geom.Domain(); amrex::GpuArray lo = box.loVect3d(); - // const auto Erads = RadSystem::ComputeThermalRadiation(T_rad_L, radBoundaries_); + // const auto Erads = RadSystem::ComputeThermalRadiation(T_rad_L, radBoundaries_); quokka::valarray const Erads = {erad_floor, EradL}; const double c_light = c; const auto Frads = Erads * c_light; @@ -154,22 +159,22 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVec // streaming inflow boundary // multigroup radiation // x1 left side boundary (Marshak) - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g]; - consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g]; - consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; - consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g]; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g]; + consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; } } // gas boundary conditions are the same everywhere const double Egas = initial_T * CV; - consVar(i, j, k, RadSystem::gasEnergy_index) = Egas; - consVar(i, j, k, RadSystem::gasDensity_index) = rho0; - consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; - consVar(i, j, k, RadSystem::x1GasMomentum_index) = 0.; - consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; - consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas; + consVar(i, j, k, RadSystem::gasDensity_index) = rho0; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; } auto problem_main() -> int @@ -187,7 +192,7 @@ auto problem_main() -> int pp.query("kappa2", kappa2); // Boundary conditions - constexpr int nvars = RadSystem::nvar_; + constexpr int nvars = RadSystem::nvar_; amrex::Vector BCs_cc(nvars); for (int n = 0; n < nvars; ++n) { BCs_cc[n].setLo(0, amrex::BCType::ext_dir); // Dirichlet x1 @@ -199,7 +204,7 @@ auto problem_main() -> int } // Problem initialization - QuokkaSimulation sim(BCs_cc); + QuokkaSimulation sim(BCs_cc); sim.radiationReconstructionOrder_ = 3; // PPM // sim.stopTime_ = tmax; // set with runtime parameters @@ -230,14 +235,14 @@ auto problem_main() -> int for (int i = 0; i < nx; ++i) { amrex::Real const x = position[i]; xs.at(i) = x; - erad1.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 0)[i]; + erad1.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 0)[i]; erad.at(i) = erad1.at(i); if (n_group_ > 1) { - erad2.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 1)[i]; + erad2.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 1)[i]; erad.at(i) += erad2.at(i); } - const double e_gas = values.at(RadSystem::gasInternalEnergy_index)[i]; - T.at(i) = quokka::EOS::ComputeTgasFromEint(rho0, e_gas); + const double e_gas = values.at(RadSystem::gasInternalEnergy_index)[i]; + T.at(i) = quokka::EOS::ComputeTgasFromEint(rho0, e_gas); T_exact.at(i) = T_end_exact; erad2_exact.at(i) = x < sim.tNew_[0] ? EradL * std::exp(-x * rho0 * kappa2) : erad_floor; diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp index 7ece497a8..3f6f475e3 100644 --- a/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp @@ -17,7 +17,7 @@ // internal headers -#include "radiation/radiation_system.hpp" +#include "radiation/radiation_dust_system.hpp" // function definitions diff --git a/src/problems/RadMarshakDustPE/CMakeLists.txt b/src/problems/RadMarshakDustPE/CMakeLists.txt new file mode 100644 index 000000000..bdcb96bc4 --- /dev/null +++ b/src/problems/RadMarshakDustPE/CMakeLists.txt @@ -0,0 +1,8 @@ +add_executable(test_radiation_marshak_dust_PE test_radiation_marshak_dust_and_PE.cpp ../../util/fextract.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radiation_marshak_dust_PE) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadiationMarshakDustPE-coupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEcoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +add_test(NAME RadiationMarshakDustPE-decoupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEdecoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp new file mode 100644 index 000000000..e86784cb4 --- /dev/null +++ b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.cpp @@ -0,0 +1,318 @@ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radiation_marshak_dust_and_PE.cpp +/// \brief Defines a test Marshak wave problem with weak coupling between dust and gas. +/// + +#include "test_radiation_marshak_dust_and_PE.hpp" +#include "AMReX.H" +#include "QuokkaSimulation.hpp" +#include "util/fextract.hpp" +#include "util/valarray.hpp" + +struct MarshakProblem { +}; + +constexpr double PE_rate = 1.0; // photoelectric heating rate in s^-1 (actual rate is PE_rate * E_FUV) +AMREX_GPU_MANAGED double kappa1 = NAN; // dust opacity at IR +AMREX_GPU_MANAGED double kappa2 = NAN; // dust opacity at FUV + +constexpr bool dust_on = 1; +constexpr bool PE_on = 1; +constexpr double gas_dust_coupling_threshold_ = 1.0e-4; + +constexpr double c = 1.0; // speed of light +constexpr double chat = 1.0; // reduced speed of light +constexpr double rho0 = 1.0; +constexpr double CV = 1.0; +constexpr double mu = 1.5 / CV; // mean molecular weight +constexpr double initial_T = 1.0; +constexpr double a_rad = 1.0; +constexpr double erad_floor = 1.0e-6; +constexpr double T_rad_L = 1.0; +constexpr double EradL = a_rad * T_rad_L * T_rad_L * T_rad_L * T_rad_L; + +constexpr int n_group_ = 2; +static constexpr amrex::GpuArray radBoundaries_{1e-10, 30, 1e4}; +static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = 1.0; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = false; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_group_; // number of radiation groups +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr int beta_order = 1; + static constexpr bool enable_dust_gas_thermal_coupling_model = dust_on; + static constexpr double energy_unit = 1.0; + static constexpr amrex::GpuArray radBoundaries = radBoundaries_; + static constexpr OpacityModel opacity_model = opacity_model_; +}; + +template <> struct ISM_Traits { + static constexpr double gas_dust_coupling_threshold = gas_dust_coupling_threshold_; + static constexpr bool enable_photoelectric_heating = PE_on; +}; + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/, + amrex::Real const /*num_density*/) -> amrex::Real +{ + // Values in cgs units from Bate & Keto (2015), Eq. 26. + // const double epsilon = 0.05; // default efficiency factor for cold molecular clouds + // const double ref_J_ISR = 5.29e-14; // reference value for the ISR in erg cm^3 + // const double coeff = 1.33e-24; + // return coeff * epsilon * num_density / ref_J_ISR; // s^-1 + + // constant rate for testing + return PE_rate; +} + +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_ + 1; ++i) { + exponents_and_values[0][i] = 0.0; + if (i == 0) { + exponents_and_values[1][i] = kappa1; + } else { + exponents_and_values[1][i] = kappa2; + } + } + return exponents_and_values; +} + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const auto Egas0 = initial_T * CV; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = erad_floor; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + 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) = Egas0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas0; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, int /*dcomp*/, int /*numcomp*/, + amrex::GeometryData const &geom, const amrex::Real /*time*/, const amrex::BCRec * /*bcr*/, + int /*bcomp*/, int /*orig_comp*/) +{ +#if (AMREX_SPACEDIM == 1) + auto i = iv.toArray()[0]; + int j = 0; + int k = 0; +#endif +#if (AMREX_SPACEDIM == 2) + auto [i, j] = iv.toArray(); + int k = 0; +#endif +#if (AMREX_SPACEDIM == 3) + auto [i, j, k] = iv.toArray(); +#endif + + amrex::Box const &box = geom.Domain(); + amrex::GpuArray lo = box.loVect3d(); + + // const auto Erads = RadSystem::ComputeThermalRadiation(T_rad_L, radBoundaries_); + quokka::valarray const Erads = {erad_floor, EradL}; + const double c_light = c; + const auto Frads = Erads * c_light; + + if (i < lo[0]) { + // streaming inflow boundary + // multigroup radiation + // x1 left side boundary (Marshak) + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g]; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g]; + consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + } + + // gas boundary conditions are the same everywhere + const double Egas = initial_T * CV; + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas; + consVar(i, j, k, RadSystem::gasDensity_index) = rho0; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; +} + +auto problem_main() -> int +{ + // Problem parameters + // const int nx = 1000; + // const double Lx = 1.0; + const double CFL_number = 0.8; + const double dt_max = 1; + const int max_timesteps = 5000; + + // read user parameters + amrex::ParmParse pp("problem"); + pp.query("kappa1", kappa1); + pp.query("kappa2", kappa2); + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + BCs_cc[n].setLo(0, amrex::BCType::ext_dir); // Dirichlet x1 + BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1 + 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 initialization + QuokkaSimulation sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + // sim.stopTime_ = tmax; // set with runtime parameters + sim.radiationCflNumber_ = CFL_number; + sim.maxDt_ = dt_max; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); + const int nx = static_cast(position.size()); + + // compute error norm + std::vector xs(nx); + std::vector T(nx); + std::vector T_exact(nx); + std::vector erad(nx); + std::vector erad1(nx); + std::vector erad2(nx); + std::vector erad1_exact(nx); + std::vector erad2_exact(nx); + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + erad1.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 0)[i]; + erad.at(i) = erad1.at(i); + if (n_group_ > 1) { + erad2.at(i) = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * 1)[i]; + erad.at(i) += erad2.at(i); + } + const double e_gas = values.at(RadSystem::gasInternalEnergy_index)[i]; + T.at(i) = quokka::EOS::ComputeTgasFromEint(rho0, e_gas); + T_exact.at(i) = x < c * sim.tNew_[0] ? initial_T + PE_rate * (sim.tNew_[0] - x / c) : initial_T; + + erad1_exact.at(i) = 0.0; + erad2_exact.at(i) = x < c * sim.tNew_[0] ? EradL : erad_floor; + } + + double err_norm = 0.; + double sol_norm = 0.; + for (int i = 1; i < nx; ++i) { // skip the first cell + err_norm += std::abs(T[i] - T_exact[i]); + err_norm += std::abs(erad1[i] - erad1_exact[i]); + err_norm += std::abs(erad2[i] - erad2_exact[i]); + sol_norm += std::abs(T_exact[i]); + sol_norm += std::abs(erad1_exact[i]); + sol_norm += std::abs(erad2_exact[i]); + } + + const double rel_err_norm = err_norm / sol_norm; + const double rel_err_tol = 0.01; + int status = 1; + if (rel_err_norm < rel_err_tol) { + status = 0; + } + amrex::Print() << "Relative L1 norm = " << rel_err_norm << std::endl; + +#ifdef HAVE_PYTHON + // Plot erad1 + matplotlibcpp::clf(); + std::map plot_args; + std::map plot_args2; + plot_args["label"] = "numerical solution"; + plot_args2["label"] = "exact solution"; + matplotlibcpp::ylim(-0.05, 1.05); + matplotlibcpp::plot(xs, erad1, plot_args); + matplotlibcpp::plot(xs, erad1_exact, plot_args2); + matplotlibcpp::xlabel("x"); + matplotlibcpp::ylabel("E_rad_group1"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("Marshak_dust test at t = {:.1f}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radiation_marshak_dust_PE_Erad1.pdf"); + + // Plot erad2 + if (n_group_ > 1) { + matplotlibcpp::clf(); + matplotlibcpp::ylim(-0.05, 1.05); + matplotlibcpp::plot(xs, erad2, plot_args); + matplotlibcpp::plot(xs, erad2_exact, plot_args2); + matplotlibcpp::xlabel("x"); + matplotlibcpp::ylabel("E_rad_group2"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("Marshak_dust test at t = {:.1f}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radiation_marshak_dust_PE_Erad2.pdf"); + } + + // plot temperature + matplotlibcpp::clf(); + matplotlibcpp::ylim(0.0, 2.1); + matplotlibcpp::plot(xs, T, plot_args); + matplotlibcpp::plot(xs, T_exact, plot_args2); + matplotlibcpp::xlabel("x"); + matplotlibcpp::ylabel("Temperature"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("Marshak_dust test at t = {:.1f}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radiation_marshak_dust_PE_temperature.pdf"); +#endif // HAVE_PYTHON + + // Cleanup and exit + amrex::Print() << "Finished." << std::endl; + return status; +} diff --git a/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.hpp b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.hpp new file mode 100644 index 000000000..dbd534a32 --- /dev/null +++ b/src/problems/RadMarshakDustPE/test_radiation_marshak_dust_and_PE.hpp @@ -0,0 +1,24 @@ +#ifndef TEST_RADIATION_MARSHAK_DUST_AND_PE_HPP_ // NOLINT +#define TEST_RADIATION_MARSHAK_DUST_AND_PE_HPP_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radiation_marshak_dust.hpp +/// \brief Defines a test Marshak wave problem with weak coupling between dust and gas. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "util/matplotlibcpp.h" +#endif +#include + +// internal headers + +#include "radiation/radiation_dust_system.hpp" + +// function definitions + +#endif // TEST_RADIATION_MARSHAK_DUST_AND_PE_HPP_ diff --git a/src/radiation/radiation_dust_system.hpp b/src/radiation/radiation_dust_system.hpp new file mode 100644 index 000000000..6d264fa82 --- /dev/null +++ b/src/radiation/radiation_dust_system.hpp @@ -0,0 +1,867 @@ +// IWYU pragma: private; include "radiation/radiation_system.hpp" +#ifndef RADIATION_DUST_SYSTEM_HPP_ +#define RADIATION_DUST_SYSTEM_HPP_ + +#define LARGE 1.e100 + +#include "radiation/radiation_system.hpp" // IWYU pragma: keep + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/, amrex::Real const num_density) + -> amrex::Real +{ + return 0.0; +} + +// Compute the Jacobian of energy update equations for the gas-dust-radiation system. The result is a struct containing the following elements: +// J00: (0, 0) component of the Jacobian matrix. = d F0 / d Egas +// F0: (0) component of the residual. = Egas residual +// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 +// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g +// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d Egas +// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g +// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual +template +AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDust( + double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, + quokka::valarray const &Src, double coeff_n, quokka::valarray const &tau, double c_v, double /*lambda_gd_time_dt*/, + quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult +{ + JacobianResult result; + + const double cscale = c_light_ / c_hat_; + + result.F0 = Egas_diff; + result.Fg = Erad_diff - (Rvec + Src); + result.Fg_abs_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] > 0.0) { + result.Fg_abs_sum += std::abs(result.Fg[g]); + result.F0 += cscale * Rvec[g]; + } + } + + // compute Jacobian elements + // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption + // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. + + auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; + + result.J00 = 1.0; + result.J0g.fillin(cscale); + const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); + dEg_dT *= d_Td_d_T; + const double dTd_dRg = -1.0 / (coeff_n * std::sqrt(T_gas)); + const auto rg = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg; + result.Jg0 = 1.0 / c_v * dEg_dT - 1.0 / cscale * rg * result.J00; + // Note that Fg is modified here, but it does not change Fg_abs_sum, which is used to check the convergence. + result.Fg = result.Fg - 1.0 / cscale * rg * result.F0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] <= 0.0) { + result.Jgg[g] = -std::numeric_limits::infinity(); + } else { + result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; + } + } + + return result; +} + +// Compute the Jacobian of energy update equations for the gas-dust-radiation system with gas and dust decoupled. The result is a struct containing the +// following elements: J00: (0, 0) component of the Jacobian matrix. = d F0 / d T_d F0: (0) component of the residual. = sum_g R_g - lambda_gd_time_dt +// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 +// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g +// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d T_d +// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g +// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual +template +AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDustDecoupled( + double /*T_gas*/, double /*T_d*/, double /*Egas_diff*/, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, + quokka::valarray const &Src, double /*coeff_n*/, quokka::valarray const &tau, double /*c_v*/, double lambda_gd_time_dt, + quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult +{ + JacobianResult result; + + result.F0 = -lambda_gd_time_dt; + result.Fg = Erad_diff - (Rvec + Src); + result.Fg_abs_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] > 0.0) { + result.F0 += Rvec[g]; + result.Fg_abs_sum += std::abs(result.Fg[g]); + } + } + + // compute Jacobian elements + // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption + // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. + + auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; + + result.J00 = 0.0; + result.J0g.fillin(1.0); + result.Jg0 = dEg_dT; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] <= 0.0) { + result.Jgg[g] = -std::numeric_limits::infinity(); + } else { + result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; + } + } + + return result; +} + +// Compute the Jacobian of energy update equations for the gas-dust-radiation system. The result is a struct containing the following elements: +// J00: (0, 0) component of the Jacobian matrix. = d F0 / d Egas +// F0: (0) component of the residual. = Egas residual +// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 +// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g +// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d Egas +// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g +// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual +template +AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDustWithPE( + double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad, quokka::valarray const &Erad0, + double PE_heating_energy_derivative, quokka::valarray const &Rvec, quokka::valarray const &Src, double coeff_n, + quokka::valarray const &tau, double c_v, double /*lambda_gd_time_dt*/, quokka::valarray const &kappaPoverE, + quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult +{ + JacobianResult result; + + const double cscale = c_light_ / c_hat_; + + result.F0 = Egas_diff; + result.Fg = Erad - Erad0 - (Rvec + Src); + result.Fg_abs_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] > 0.0) { + result.Fg_abs_sum += std::abs(result.Fg[g]); + result.F0 += cscale * Rvec[g]; + } + } + result.F0 -= PE_heating_energy_derivative * Erad[nGroups_ - 1]; + + // const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivativeMultiGroup(T_d, radBoundaries_g_copy); + AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); + + // compute Jacobian elements + // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption + // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. + + const auto d_Eg_d_Rg = -1.0 * kappaPoverE / tau; + + result.J00 = 1.0; + result.J0g.fillin(cscale); + if (tau[nGroups_ - 1] <= 0.0) { + result.J0g[nGroups_ - 1] = LARGE; + } else { + result.J0g[nGroups_ - 1] -= PE_heating_energy_derivative * d_Eg_d_Rg[nGroups_ - 1]; + } + const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); + const auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t * d_Td_d_T; + const double dTd_dRg = -1.0 / (coeff_n * std::sqrt(T_gas)); + const auto rg = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg; + result.Jg0 = 1.0 / c_v * dEg_dT - 1.0 / cscale * rg * result.J00; + // Note that Fg is modified here, but it does not change Fg_abs_sum, which is used to check the convergence. + result.Fg = result.Fg - 1.0 / cscale * rg * result.F0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] <= 0.0) { + result.Jgg[g] = -LARGE; + } else { + result.Jgg[g] = d_Eg_d_Rg[g] - 1.0; + } + } + result.Jgg[nGroups_ - 1] += rg[nGroups_ - 1] / cscale * PE_heating_energy_derivative * d_Eg_d_Rg[nGroups_ - 1]; + result.Jg1 = rg - 1.0 / cscale * rg * result.J0g[nGroups_ - 1]; // note that this is the (nGroups_ - 1)th column, except for the (nGroups_ - 1)th row + + return result; +} + +// 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] +// [ai0 aii] [xi] [yi] +// for x0 and xi, where a0i = (a01, a02, a03, ...); ai0 = (a10, a20, a30, ...); aii = (a11, a22, a33, ...), xi = (x1, x2, x3, ...), yi = (y1, y2, y3, ...) +template +AMREX_GPU_HOST_DEVICE void RadSystem::SolveLinearEqsWithLastColumn(JacobianResult const &jacobian, double &x0, + quokka::valarray &xi) +{ + // Note that the following routine only works when the FUV group is the last group, i.e., pe_index = nGroups_ - 1 + + // note that jacobian.Jgg[pe_index] is the right value for a[pe_index][pe_index], while jacobian.Jg1[pe_index] is NOT. + const int pe_index = nGroups_ - 1; + const auto ratios = jacobian.J0g / jacobian.Jgg; + + const auto a00_new = jacobian.J00 - sum(ratios * jacobian.Jg0); + const auto y0_new = jacobian.F0 - sum(ratios * jacobian.Fg); + auto a01_new = jacobian.J0g[pe_index] - sum(ratios * jacobian.Jg1); + // re-accounting for the pe_index'th element of jacobian.Jg1 + a01_new = a01_new + ratios[pe_index] * jacobian.Jg1[pe_index] - ratios[pe_index] * jacobian.Jgg[pe_index]; + const auto a10 = jacobian.Jg0[pe_index]; + const auto a11 = jacobian.Jgg[pe_index]; + const auto y1 = jacobian.Fg[pe_index]; + // solve linear equations [[a00_new, a01_new], [a10, a11]] [[x0], [xi[pe_index]]] = [y0_new, y1] + x0 = (y0_new - a01_new / a11 * y1) / (a00_new - a01_new / a11 * a10); + const auto x1 = (y1 - a10 * x0) / a11; + xi[pe_index] = x1; + // xi = (jacobian.Fg - jacobian.Jg0 * x0) / jacobian.Jgg; + for (int g = 0; g < pe_index; ++g) { + xi[g] = (jacobian.Fg[g] - jacobian.Jg0[g] * x0 - jacobian.Jg1[g] * x1) / jacobian.Jgg[g]; + } + x0 *= -1.0; + xi *= -1.0; +} + +template +AMREX_GPU_DEVICE auto RadSystem::SolveGasDustRadiationEnergyExchange( + double const Egas0, quokka::valarray const &Erad0Vec, double const rho, double const coeff_n, double const dt, + amrex::GpuArray const &massScalars, int const n_outer_iter, quokka::valarray const &work, + quokka::valarray const &vel_times_F, quokka::valarray const &Src, + amrex::GpuArray const &rad_boundaries, int *p_iteration_counter, int *p_iteration_failure_counter) -> NewtonIterationResult +{ + // 1. Compute energy exchange + + // BEGIN NEWTON-RAPHSON LOOP + // Define the source term: S = dt chat gamma rho (kappa_P B - kappa_E E) + dt chat c^-2 gamma rho kappa_F v * F_i, where gamma = + // 1 / sqrt(1 - v^2 / c^2) is the Lorentz factor. Solve for the new radiation energy and gas internal energy using a + // Newton-Raphson method using the base variables (Egas, D_0, D_1, + // ...), where D_i = R_i / tau_i^(t) and tau_i^(t) = dt * chat * gamma * rho * kappa_{P,i}^(t) is the optical depth across chat + // * dt for group i at time t. Compared with the old base (Egas, Erad_0, Erad_1, ...), this new base is more stable and + // converges faster. Furthermore, the PlanckOpacityTempDerivative term is not needed anymore since we assume d/dT (kappa_P / + // kappa_E) = 0 in the calculation of the Jacobian. Note that this assumption only affects the convergence rate of the + // Newton-Raphson iteration and does not affect the result at all once the iteration is converged. + // + // The Jacobian of F(E_g, D_i) is + // + // dF_G / dE_g = 1 + // dF_G / dD_i = c / chat * tau0_i + // dF_{D,i} / dE_g = 1 / (chat * C_v) * (kappa_{P,i} / kappa_{E,i}) * d/dT (4 \pi B_i) + // dF_{D,i} / dD_i = - (1 / (chat * dt * rho * kappa_{E,i}) + 1) * tau0_i = - ((1 / tau_i)(kappa_Pi / kappa_Ei) + 1) * tau0_i + + const double c = c_light_; // make a copy of c_light_ to avoid compiler error "undefined in device code" + const double chat = c_hat_; + const double cscale = c / chat; + + int dust_model = 1; + double T_d0 = NAN; + double lambda_gd_times_dt = NAN; + const double T_gas0 = quokka::EOS::ComputeTgasFromEint(rho, Egas0, massScalars); + AMREX_ASSERT(T_gas0 >= 0.); + T_d0 = ComputeDustTemperatureBateKeto(T_gas0, T_gas0, rho, Erad0Vec, coeff_n, dt, NAN, 0, rad_boundaries); + AMREX_ASSERT_WITH_MESSAGE(T_d0 >= 0., "Dust temperature is negative!"); + if (T_d0 < 0.0) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT + } + + const double max_Gamma_gd = coeff_n * std::max(std::sqrt(T_gas0) * T_gas0, std::sqrt(T_d0) * T_d0); + if (cscale * max_Gamma_gd < ISM_Traits::gas_dust_coupling_threshold * Egas0) { + dust_model = 2; + lambda_gd_times_dt = coeff_n * std::sqrt(T_gas0) * (T_gas0 - T_d0); + } + + // const double Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); + double Etot0 = NAN; + if (dust_model == 1) { + Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); + } else { + // for dust_model == 2 (decoupled gas and dust), Egas0 is not involved in the iteration + const double fourPiBoverC = sum(ComputeThermalRadiationMultiGroup(T_d0, rad_boundaries)); + Etot0 = std::abs(lambda_gd_times_dt) + fourPiBoverC + (sum(Erad0Vec) + sum(Src)); + } + + double T_gas = NAN; + double T_d = NAN; + double delta_x = NAN; + quokka::valarray delta_R{}; + quokka::valarray F_D{}; + quokka::valarray Rvec{}; + quokka::valarray tau0{}; // optical depth across c * dt at old state + quokka::valarray tau{}; // optical depth across c * dt at new state + quokka::valarray work_local{}; // work term used in the Newton-Raphson iteration of the current outer iteration + quokka::valarray fourPiBoverC{}; + amrex::GpuArray rad_boundary_ratios{}; + amrex::GpuArray, 2> kappa_expo_and_lower_value{}; + OpacityTerms opacity_terms{}; + + // fill kappa_expo_and_lower_value with NAN to get warned when there are uninitialized values + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < nGroups_ + 1; ++j) { + kappa_expo_and_lower_value[i][j] = NAN; + } + } + + if constexpr (!(opacity_model_ == OpacityModel::piecewise_constant_opacity)) { + for (int g = 0; g < nGroups_; ++g) { + rad_boundary_ratios[g] = rad_boundaries[g + 1] / rad_boundaries[g]; + } + } + + // define a list of alpha_quant for the model PPL_opacity_fixed_slope_spectrum + amrex::GpuArray alpha_quant_minus_one{}; + if constexpr ((opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) || + (gamma_ == 1.0 && opacity_model_ == OpacityModel::PPL_opacity_full_spectrum)) { + if constexpr (!special_edge_bin_slopes) { + for (int g = 0; g < nGroups_; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } else { + alpha_quant_minus_one[0] = 2.0; + alpha_quant_minus_one[nGroups_ - 1] = -4.0; + for (int g = 1; g < nGroups_ - 1; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } + } + + double Egas_guess = Egas0; + auto EradVec_guess = Erad0Vec; + + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + + if (dust_model == 2) { + Egas_guess = Egas0 - cscale * lambda_gd_times_dt; // update Egas_guess once for all + } + + const double resid_tol = 1.0e-11; // 1.0e-15; + const int maxIter = 100; + int n = 0; + for (; n < maxIter; ++n) { + // 1. Compute dust temperature + // If the dust model is turned off, ComputeDustTemperature should be a function that returns T_gas. + + if (n > 0) { + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + } + + if (dust_model == 1) { + if (n == 0) { + T_d = T_d0; + } else { + T_d = T_gas - sum(Rvec) / (coeff_n * std::sqrt(T_gas)); + } + } else { + if (n == 0) { + T_d = T_d0; + } + } + AMREX_ASSERT_WITH_MESSAGE(T_d >= 0., "Dust temperature is negative! Consider increasing ISM_Traits::gas_dust_coupling_threshold"); + if (T_d < 0.0) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT + } + + // 2. Compute kappaP and kappaE at dust temperature + + fourPiBoverC = ComputeThermalRadiationMultiGroup(T_d, rad_boundaries); + + opacity_terms = ComputeModelDependentKappaEAndKappaP(T_d, rho, rad_boundaries, rad_boundary_ratios, fourPiBoverC, EradVec_guess, n, + opacity_terms.alpha_E, opacity_terms.alpha_P); + + if (n == 0) { + // Compute kappaF and the delta_nu_kappa_B term. kappaF is used to compute the work term. + // Will update opacity_terms in place + ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place + } + + // 3. In the first loop, calculate kappaF, work, tau0, R + + if (n == 0) { + + if constexpr ((beta_order_ == 1) && (include_work_term_in_source)) { + // compute the work term at the old state + // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); + if (n_outer_iter == 0) { + for (int g = 0; g < nGroups_; ++g) { + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + work_local[g] = vel_times_F[g] * opacity_terms.kappaF[g] * chat / (c * c) * dt; + } else { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(rad_boundaries, rho, T_d); + work_local[g] = vel_times_F[g] * opacity_terms.kappaF[g] * chat / (c * c) * dt * + (1.0 + kappa_expo_and_lower_value[0][g]); + } + } + } else { + // If n_outer_iter > 0, use the work term from the previous outer iteration, which is passed as the parameter 'work' + work_local = work; + } + } else { + work_local.fillin(0.0); + } + + tau0 = dt * rho * opacity_terms.kappaP * chat; + tau = tau0; + Rvec = (fourPiBoverC - EradVec_guess / opacity_terms.kappaPoverE) * tau0 + work_local; + if constexpr (use_D_as_base) { + // tau0 is used as a scaling factor for Rvec + for (int g = 0; g < nGroups_; ++g) { + if (tau0[g] <= 1.0) { + tau0[g] = 1.0; + } + } + } + } else { // in the second and later loops, calculate tau and E (given R) + tau = dt * rho * opacity_terms.kappaP * chat; + for (int g = 0; g < nGroups_; ++g) { + // If tau = 0.0, Erad_guess shouldn't change + if (tau[g] > 0.0) { + EradVec_guess[g] = opacity_terms.kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work_local[g]) / tau[g]); + if constexpr (force_rad_floor_in_iteration) { + if (EradVec_guess[g] < 0.0) { + Egas_guess -= cscale * (Erad_floor_ - EradVec_guess[g]); + EradVec_guess[g] = Erad_floor_; + } + } + } + } + } + + const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivativeMultiGroup(T_d, rad_boundaries); + AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); + const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T + + const auto Egas_diff = Egas_guess - Egas0; + const auto Erad_diff = EradVec_guess - Erad0Vec; + + JacobianResult jacobian; + + if (dust_model == 1) { + jacobian = ComputeJacobianForGasAndDust(T_gas, T_d, Egas_diff, Erad_diff, Rvec, Src, coeff_n, tau, c_v, lambda_gd_times_dt, + opacity_terms.kappaPoverE, d_fourpiboverc_d_t); + } else { + jacobian = ComputeJacobianForGasAndDustDecoupled(T_gas, T_d, Egas_diff, Erad_diff, Rvec, Src, coeff_n, tau, c_v, lambda_gd_times_dt, + opacity_terms.kappaPoverE, d_fourpiboverc_d_t); + } + + if constexpr (use_D_as_base) { + jacobian.J0g = jacobian.J0g * tau0; + jacobian.Jgg = jacobian.Jgg * tau0; + } + + // check relative convergence of the residuals + if ((std::abs(jacobian.F0 / Etot0) < resid_tol) && (cscale * jacobian.Fg_abs_sum / Etot0 < resid_tol)) { + break; + } + +#if 0 + // For debugging: print (Egas0, Erad0Vec, tau0), which defines the initial condition for a Newton-Raphson iteration + if (n == 0) { + std::cout << "Egas0 = " << Egas0 << ", Erad0Vec = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << Erad0Vec[g] << ", "; + } + std::cout << "], tau0 = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << tau0[g] << ", "; + } + std::cout << "]"; + std::cout << "; C_V = " << c_v << ", a_rad = " << radiation_constant_ << ", coeff_n = " << coeff_n << "\n"; + } else if (n >= 0) { + std::cout << "n = " << n << ", Egas_guess = " << Egas_guess << ", EradVec_guess = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << EradVec_guess[g] << ", "; + } + std::cout << "], tau = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << tau[g] << ", "; + } + std::cout << "]"; + std::cout << ", F_G = " << jacobian.F0 << ", F_D_abs_sum = " << jacobian.Fg_abs_sum << ", Etot0 = " << Etot0 << "\n"; + } +#endif + + // update variables + RadSystem::SolveLinearEqs(jacobian, delta_x, delta_R); // This is modify delta_x and delta_R in place + AMREX_ASSERT(!std::isnan(delta_x)); + AMREX_ASSERT(!delta_R.hasnan()); + + // Update independent variables (Egas_guess, Rvec) + // enable_dE_constrain is used to prevent the gas temperature from dropping/increasing below/above the radiation + // temperature + if (dust_model == 2) { + T_d += delta_x; + Rvec += delta_R; + } else { + const double T_rad = std::sqrt(std::sqrt(sum(EradVec_guess) / radiation_constant_)); + if (enable_dE_constrain && delta_x / c_v > std::max(T_gas, T_rad)) { + Egas_guess = quokka::EOS::ComputeEintFromTgas(rho, T_rad); + // Rvec.fillin(0.0); + } else { + Egas_guess += delta_x; + if constexpr (use_D_as_base) { + Rvec += tau0 * delta_R; + } else { + Rvec += delta_R; + } + } + } + + // check relative and absolute convergence of E_r + // if (std::abs(deltaEgas / Egas_guess) < 1e-7) { + // break; + // } + } // END NEWTON-RAPHSON LOOP + + AMREX_ASSERT(Egas_guess > 0.0); + AMREX_ASSERT(min(EradVec_guess) >= 0.0); + + AMREX_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration for matter-radiation coupling failed to converge!"); + if (n >= maxIter) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[0], 1); // NOLINT + } + + amrex::Gpu::Atomic::Add(&p_iteration_counter[0], 1); // total number of radiation updates. NOLINT + amrex::Gpu::Atomic::Add(&p_iteration_counter[1], n + 1); // total number of Newton-Raphson iterations. NOLINT + amrex::Gpu::Atomic::Max(&p_iteration_counter[2], n + 1); // maximum number of Newton-Raphson iterations. NOLINT + if (dust_model == 2) { + amrex::Gpu::Atomic::Add(&p_iteration_counter[3], 1); // total number of decoupled gas-dust iterations. NOLINT + } + + NewtonIterationResult result; + + if (n > 0) { + // calculate kappaF since the temperature has changed + // Will update opacity_terms in place + ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place + } + + result.Egas = Egas_guess; + result.EradVec = EradVec_guess; + result.work = work_local; + result.T_gas = T_gas; + result.T_d = T_d; + result.opacity_terms = opacity_terms; + return result; +} + +template +AMREX_GPU_DEVICE auto RadSystem::SolveGasDustRadiationEnergyExchangeWithPE( + double const Egas0, quokka::valarray const &Erad0Vec, double const rho, double const coeff_n, double const dt, + amrex::GpuArray const &massScalars, int const n_outer_iter, quokka::valarray const &work, + quokka::valarray const &vel_times_F, quokka::valarray const &Src, + amrex::GpuArray const &rad_boundaries, int *p_iteration_counter, int *p_iteration_failure_counter) -> NewtonIterationResult +{ + // 1. Compute energy exchange + + // BEGIN NEWTON-RAPHSON LOOP + // Define the source term: S = dt chat gamma rho (kappa_P B - kappa_E E) + dt chat c^-2 gamma rho kappa_F v * F_i, where gamma = + // 1 / sqrt(1 - v^2 / c^2) is the Lorentz factor. Solve for the new radiation energy and gas internal energy using a + // Newton-Raphson method using the base variables (Egas, D_0, D_1, + // ...), where D_i = R_i / tau_i^(t) and tau_i^(t) = dt * chat * gamma * rho * kappa_{P,i}^(t) is the optical depth across chat + // * dt for group i at time t. Compared with the old base (Egas, Erad_0, Erad_1, ...), this new base is more stable and + // converges faster. Furthermore, the PlanckOpacityTempDerivative term is not needed anymore since we assume d/dT (kappa_P / + // kappa_E) = 0 in the calculation of the Jacobian. Note that this assumption only affects the convergence rate of the + // Newton-Raphson iteration and does not affect the result at all once the iteration is converged. + // + // The Jacobian of F(E_g, D_i) is + // + // dF_G / dE_g = 1 + // dF_G / dD_i = c / chat * tau0_i + // dF_{D,i} / dE_g = 1 / (chat * C_v) * (kappa_{P,i} / kappa_{E,i}) * d/dT (4 \pi B_i) + // dF_{D,i} / dD_i = - (1 / (chat * dt * rho * kappa_{E,i}) + 1) * tau0_i = - ((1 / tau_i)(kappa_Pi / kappa_Ei) + 1) * tau0_i + + const double c = c_light_; // make a copy of c_light_ to avoid compiler error "undefined in device code" + const double chat = c_hat_; + const double cscale = c / chat; + + int dust_model = 1; + double T_d0 = NAN; + double lambda_gd_times_dt = NAN; + const double T_gas0 = quokka::EOS::ComputeTgasFromEint(rho, Egas0, massScalars); + AMREX_ASSERT(T_gas0 >= 0.); + T_d0 = ComputeDustTemperatureBateKeto(T_gas0, T_gas0, rho, Erad0Vec, coeff_n, dt, NAN, 0, rad_boundaries); + AMREX_ASSERT_WITH_MESSAGE(T_d0 >= 0., "Dust temperature is negative!"); + if (T_d0 < 0.0) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT + } + + const double max_Gamma_gd = coeff_n * std::max(std::sqrt(T_gas0) * T_gas0, std::sqrt(T_d0) * T_d0); + if (cscale * max_Gamma_gd < ISM_Traits::gas_dust_coupling_threshold * Egas0) { + dust_model = 2; + lambda_gd_times_dt = coeff_n * std::sqrt(T_gas0) * (T_gas0 - T_d0); + } + + // const double Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); + double Etot0 = NAN; + if (dust_model == 1) { + Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); + } else { + // for dust_model == 2 (decoupled gas and dust), Egas0 is not involved in the iteration + Etot0 = std::abs(lambda_gd_times_dt) + (sum(Erad0Vec) + sum(Src)); + } + + double T_gas = NAN; + double T_d = NAN; + double delta_x = NAN; + quokka::valarray delta_R{}; + quokka::valarray F_D{}; + quokka::valarray Rvec{}; + quokka::valarray tau0{}; // optical depth across c * dt at old state + quokka::valarray tau{}; // optical depth across c * dt at new state + quokka::valarray work_local{}; // work term used in the Newton-Raphson iteration of the current outer iteration + quokka::valarray fourPiBoverC{}; + amrex::GpuArray rad_boundary_ratios{}; + amrex::GpuArray, 2> kappa_expo_and_lower_value{}; + OpacityTerms opacity_terms{}; + + // fill kappa_expo_and_lower_value with NAN to avoid mistakes caused by uninitialized values + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < nGroups_ + 1; ++j) { + kappa_expo_and_lower_value[i][j] = NAN; + } + } + + if constexpr (!(opacity_model_ == OpacityModel::piecewise_constant_opacity)) { + for (int g = 0; g < nGroups_; ++g) { + rad_boundary_ratios[g] = rad_boundaries[g + 1] / rad_boundaries[g]; + } + } + + // define a list of alpha_quant for the model PPL_opacity_fixed_slope_spectrum + amrex::GpuArray alpha_quant_minus_one{}; + if constexpr ((opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) || + (gamma_ == 1.0 && opacity_model_ == OpacityModel::PPL_opacity_full_spectrum)) { + if constexpr (!special_edge_bin_slopes) { + for (int g = 0; g < nGroups_; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } else { + alpha_quant_minus_one[0] = 2.0; + alpha_quant_minus_one[nGroups_ - 1] = -4.0; + for (int g = 1; g < nGroups_ - 1; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } + } + + double Egas_guess = Egas0; + auto EradVec_guess = Erad0Vec; + + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + + if (dust_model == 2) { + Egas_guess = Egas0 - cscale * lambda_gd_times_dt; // update Egas_guess once for all + } + + // phtoelectric heating + const double num_den = rho / mean_molecular_mass_; + const double PE_heating_energy_derivative = dt * DefinePhotoelectricHeatingE1Derivative(T_gas, num_den); + + const double resid_tol = 1.0e-11; // 1.0e-15; + const int maxIter = 100; + int n = 0; + for (; n < maxIter; ++n) { + // 1. Compute dust temperature + // If the dust model is turned off, ComputeDustTemperature should be a function that returns T_gas. + + if (n > 0) { + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + } + + if (dust_model == 1) { + if (n == 0) { + T_d = T_d0; + } else { + T_d = T_gas - sum(Rvec) / (coeff_n * std::sqrt(T_gas)); + } + } else { + if (n == 0) { + T_d = T_d0; + } + } + AMREX_ASSERT_WITH_MESSAGE(T_d >= 0., "Dust temperature is negative! Consider increasing ISM_Traits::gas_dust_coupling_threshold"); + if (T_d < 0.0) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT + } + + // 2. Compute kappaP and kappaE at dust temperature + + fourPiBoverC = ComputeThermalRadiationMultiGroup(T_d, rad_boundaries); + + opacity_terms = ComputeModelDependentKappaEAndKappaP(T_d, rho, rad_boundaries, rad_boundary_ratios, fourPiBoverC, EradVec_guess, n, + opacity_terms.alpha_E, opacity_terms.alpha_P); + + if (n == 0) { + // Compute kappaF and the delta_nu_kappa_B term. kappaF is used to compute the work term. + // Only the last two arguments (kappaFVec, delta_nu_kappa_B_at_edge) are modified in this function. + ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place + } + + // 3. In the first loop, calculate kappaF, work, tau0, R + + if (n == 0) { + + if constexpr ((beta_order_ == 1) && (include_work_term_in_source)) { + // compute the work term at the old state + // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); + if (n_outer_iter == 0) { + for (int g = 0; g < nGroups_; ++g) { + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + work_local[g] = vel_times_F[g] * opacity_terms.kappaF[g] * chat / (c * c) * dt; + } else { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(rad_boundaries, rho, T_d); + work_local[g] = vel_times_F[g] * opacity_terms.kappaF[g] * chat / (c * c) * dt * + (1.0 + kappa_expo_and_lower_value[0][g]); + } + } + } else { + // If n_outer_iter > 0, use the work term from the previous outer iteration, which is passed as the parameter 'work' + work_local = work; + } + } else { + work_local.fillin(0.0); + } + + tau0 = dt * rho * opacity_terms.kappaP * chat; + tau = tau0; + Rvec = (fourPiBoverC - EradVec_guess / opacity_terms.kappaPoverE) * tau0 + work_local; + if constexpr (use_D_as_base) { + // tau0 is used as a scaling factor for Rvec + for (int g = 0; g < nGroups_; ++g) { + if (tau0[g] <= 1.0) { + tau0[g] = 1.0; + } + } + } + } else { // in the second and later loops, calculate tau and E (given R) + tau = dt * rho * opacity_terms.kappaP * chat; + for (int g = 0; g < nGroups_; ++g) { + // If tau = 0.0, Erad_guess shouldn't change + if (tau[g] > 0.0) { + EradVec_guess[g] = opacity_terms.kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work_local[g]) / tau[g]); + if constexpr (force_rad_floor_in_iteration) { + if (EradVec_guess[g] < 0.0) { + Egas_guess -= cscale * (Erad_floor_ - EradVec_guess[g]); + EradVec_guess[g] = Erad_floor_; + } + } + } + } + } + + const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivativeMultiGroup(T_d, rad_boundaries); + AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); + const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T + + const auto Egas_diff = Egas_guess - Egas0; + const auto Erad_diff = EradVec_guess - Erad0Vec; + + JacobianResult jacobian; + + if (dust_model == 1) { + jacobian = ComputeJacobianForGasAndDustWithPE(T_gas, T_d, Egas_diff, EradVec_guess, Erad0Vec, PE_heating_energy_derivative, Rvec, Src, + coeff_n, tau, c_v, lambda_gd_times_dt, opacity_terms.kappaPoverE, d_fourpiboverc_d_t); + } else { + jacobian = ComputeJacobianForGasAndDustDecoupled(T_gas, T_d, Egas_diff, Erad_diff, Rvec, Src, coeff_n, tau, c_v, lambda_gd_times_dt, + opacity_terms.kappaPoverE, d_fourpiboverc_d_t); + } + + if constexpr (use_D_as_base) { + jacobian.J0g = jacobian.J0g * tau0; + jacobian.Jgg = jacobian.Jgg * tau0; + } + + // check relative convergence of the residuals + if ((std::abs(jacobian.F0 / Etot0) < resid_tol) && (cscale * jacobian.Fg_abs_sum / Etot0 < resid_tol)) { + break; + } + +#if 0 + // For debugging: print (Egas0, Erad0Vec, tau0), which defines the initial condition for a Newton-Raphson iteration + if (n == 0) { + std::cout << "Egas0 = " << Egas0 << ", Erad0Vec = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << Erad0Vec[g] << ", "; + } + std::cout << "], tau0 = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << tau0[g] << ", "; + } + std::cout << "]"; + std::cout << "; C_V = " << c_v << ", a_rad = " << radiation_constant_ << ", coeff_n = " << coeff_n << "\n"; + } else if (n >= 0) { + std::cout << "n = " << n << ", Egas_guess = " << Egas_guess << ", EradVec_guess = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << EradVec_guess[g] << ", "; + } + std::cout << "], tau = ["; + for (int g = 0; g < nGroups_; ++g) { + std::cout << tau[g] << ", "; + } + std::cout << "]"; + std::cout << ", F_G = " << jacobian.F0 << ", F_D_abs_sum = " << jacobian.Fg_abs_sum << ", Etot0 = " << Etot0 << "\n"; + } +#endif + + // update variables + RadSystem::SolveLinearEqsWithLastColumn(jacobian, delta_x, delta_R); // This is modify delta_x and delta_R in place + AMREX_ASSERT(!std::isnan(delta_x)); + AMREX_ASSERT(!delta_R.hasnan()); + + // Update independent variables (Egas_guess, Rvec) + // enable_dE_constrain is used to prevent the gas temperature from dropping/increasing below/above the radiation + // temperature + if (dust_model == 2) { + T_d += delta_x; + Rvec += delta_R; + } else { + const double T_rad = std::sqrt(std::sqrt(sum(EradVec_guess) / radiation_constant_)); + if (enable_dE_constrain && delta_x / c_v > std::max(T_gas, T_rad)) { + Egas_guess = quokka::EOS::ComputeEintFromTgas(rho, T_rad); + // Rvec.fillin(0.0); + } else { + Egas_guess += delta_x; + if constexpr (use_D_as_base) { + Rvec += tau0 * delta_R; + } else { + Rvec += delta_R; + } + } + } + + // check relative and absolute convergence of E_r + // if (std::abs(deltaEgas / Egas_guess) < 1e-7) { + // break; + // } + } // END NEWTON-RAPHSON LOOP + + if (dust_model == 2) { + Egas_guess += PE_heating_energy_derivative * EradVec_guess[nGroups_ - 1]; + } + + AMREX_ASSERT(Egas_guess > 0.0); + AMREX_ASSERT(min(EradVec_guess) >= 0.0); + + AMREX_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); + if (n >= maxIter) { + amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[0], 1); // NOLINT + } + + amrex::Gpu::Atomic::Add(&p_iteration_counter[0], 1); // total number of radiation updates. NOLINT + amrex::Gpu::Atomic::Add(&p_iteration_counter[1], n + 1); // total number of Newton-Raphson iterations. NOLINT + amrex::Gpu::Atomic::Max(&p_iteration_counter[2], n + 1); // maximum number of Newton-Raphson iterations. NOLINT + if (dust_model == 2) { + amrex::Gpu::Atomic::Add(&p_iteration_counter[3], 1); // total number of decoupled gas-dust iterations. NOLINT + } + + NewtonIterationResult result; + + if (n > 0) { + // calculate kappaF since the temperature has changed + // Will update opacity_terms in place + ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place + } + + result.Egas = Egas_guess; + result.EradVec = EradVec_guess; + result.work = work_local; + result.T_gas = T_gas; + result.T_d = T_d; + result.opacity_terms = opacity_terms; + return result; +} + +#endif // RADIATION_DUST_SYSTEM_HPP_ \ No newline at end of file diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 36b5c2d64..678cd0232 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -84,6 +84,7 @@ template struct RadSystem_Traits { // template struct ISM_Traits { static constexpr double gas_dust_coupling_threshold = 1.0e-6; + static constexpr bool enable_photoelectric_heating = false; }; // A struct to hold the results of the ComputeRadPressure function. @@ -124,6 +125,7 @@ template struct JacobianResult { quokka::valarray::nGroups> J0g; // (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups quokka::valarray::nGroups> Jg0; // (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups quokka::valarray::nGroups> Jgg; // (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups + quokka::valarray::nGroups> Jg1; // (g, 1) components of the Jacobian matrix, g = 1, 2, ..., nGroups quokka::valarray::nGroups> Fg; // (g) components of the residual, g = 1, 2, ..., nGroups }; @@ -192,6 +194,7 @@ template class RadSystem : public HyperbolicSystem::beta_order; static constexpr bool enable_dust_gas_thermal_coupling_model_ = RadSystem_Traits::enable_dust_gas_thermal_coupling_model; + static constexpr bool enable_photoelectric_heating_ = ISM_Traits::enable_photoelectric_heating; static constexpr int nGroups_ = Physics_Traits::nGroups; static constexpr amrex::GpuArray radBoundaries_ = []() constexpr { @@ -300,6 +303,9 @@ template class RadSystem : public HyperbolicSystem const &jacobian, double &x0, quokka::valarray &xi); + AMREX_GPU_HOST_DEVICE static void SolveLinearEqsWithLastColumn(JacobianResult const &jacobian, 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) -> std::tuple; @@ -328,9 +334,11 @@ template class RadSystem : public HyperbolicSystem const &rad_boundaries = amrex::GpuArray{}, amrex::GpuArray const &rad_boundary_ratios = amrex::GpuArray{}) -> double; - AMREX_GPU_DEVICE static auto ComputeJacobianForGas(double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, + AMREX_GPU_HOST_DEVICE static auto DefinePhotoelectricHeatingE1Derivative(amrex::Real temperature, amrex::Real num_density) -> amrex::Real; + + AMREX_GPU_DEVICE static auto ComputeJacobianForGas(double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, quokka::valarray const &Src, - double coeff_n, quokka::valarray const &tau, double c_v, double lambda_gd_time_dt, + quokka::valarray const &tau, double c_v, quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult; @@ -357,14 +365,38 @@ template class RadSystem : public HyperbolicSystem const &alpha_E, amrex::GpuArray const &alpha_P) -> OpacityTerms; - template + AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDustWithPE( + double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad, quokka::valarray const &Erad0, + double PE_heating_energy_derivative, quokka::valarray const &Rvec, quokka::valarray const &Src, double coeff_n, + quokka::valarray const &tau, double c_v, double lambda_gd_time_dt, quokka::valarray const &kappaPoverE, + quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult; + + AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDustDecoupledWithPE( + double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, + quokka::valarray const &Src, double coeff_n, quokka::valarray const &tau, double c_v, double lambda_gd_time_dt, + quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult; + + AMREX_GPU_DEVICE static auto + SolveGasRadiationEnergyExchange(double Egas0, quokka::valarray const &Erad0Vec, double rho, double dt, + amrex::GpuArray const &massScalars, int n_outer_iter, quokka::valarray const &work, + quokka::valarray const &vel_times_F, quokka::valarray const &Src, + amrex::GpuArray const &rad_boundaries, int *p_iteration_counter, int *p_iteration_failure_counter) + -> NewtonIterationResult; + + AMREX_GPU_DEVICE static auto SolveGasDustRadiationEnergyExchange(double Egas0, quokka::valarray const &Erad0Vec, double rho, + double coeff_n, double dt, amrex::GpuArray const &massScalars, + int n_outer_iter, quokka::valarray const &work, + quokka::valarray const &vel_times_F, + quokka::valarray const &Src, + amrex::GpuArray const &rad_boundaries, int *p_iteration_counter, + int *p_iteration_failure_counter) -> NewtonIterationResult; + AMREX_GPU_DEVICE static auto - SolveMatterRadiationEnergyExchange(double Egas0, quokka::valarray const &Erad0Vec, double rho, double T_d0, int dust_model, - double coeff_n, double lambda_gd_times_dt, double dt, amrex::GpuArray const &massScalars, - int n_outer_iter, quokka::valarray const &work, - quokka::valarray const &vel_times_F, quokka::valarray const &Src, - amrex::GpuArray const &rad_boundaries, JacobianFunc ComputeJacobian, int *p_iteration_counter, - int *p_iteration_failure_counter) -> NewtonIterationResult; + SolveGasDustRadiationEnergyExchangeWithPE(double Egas0, quokka::valarray const &Erad0Vec, double rho, double coeff_n, double dt, + amrex::GpuArray const &massScalars, int n_outer_iter, + quokka::valarray const &work, quokka::valarray const &vel_times_F, + quokka::valarray const &Src, amrex::GpuArray const &rad_boundaries, + int *p_iteration_counter, int *p_iteration_failure_counter) -> NewtonIterationResult; template AMREX_GPU_DEVICE static auto ComputeCellOpticalDepth(const quokka::Array4View &consVar, @@ -1367,7 +1399,7 @@ AMREX_GPU_DEVICE auto RadSystem::ComputeDustTemperatureBateKeto(doubl const double LHS = c_hat_ * dt * rho * sum(kappaEVec * Erad - kappaPVec * fourPiBoverC) + N_d * std::sqrt(T_gas) * (T_gas - T_d); - if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { // TODO: remove abs here + if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { break; } diff --git a/src/radiation/source_terms_multi_group.hpp b/src/radiation/source_terms_multi_group.hpp index 56a329ff3..f8069442e 100644 --- a/src/radiation/source_terms_multi_group.hpp +++ b/src/radiation/source_terms_multi_group.hpp @@ -4,162 +4,6 @@ #include "radiation/radiation_system.hpp" // IWYU pragma: keep -// Compute the Jacobian of energy update equations for the gas-radiation system. The result is a struct containing the following elements: -// J00: (0, 0) component of the Jacobian matrix. = d F0 / d Egas -// F0: (0) component of the residual. = Egas residual -// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 -// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g -// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d Egas -// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g -// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual -template -AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGas(double /*T_gas*/, double /*T_d*/, double Egas_diff, - quokka::valarray const &Erad_diff, - quokka::valarray const &Rvec, quokka::valarray const &Src, - double /*coeff_n*/, quokka::valarray const &tau, double c_v, - double /*lambda_gd_time_dt*/, quokka::valarray const &kappaPoverE, - quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult -{ - JacobianResult result; - - const double cscale = c_light_ / c_hat_; - - result.F0 = Egas_diff; - result.Fg = Erad_diff - (Rvec + Src); - result.Fg_abs_sum = 0.0; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] > 0.0) { - result.Fg_abs_sum += std::abs(result.Fg[g]); - result.F0 += cscale * Rvec[g]; - } - } - - // const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivativeMultiGroup(T_d, radBoundaries_g_copy); - AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); - - // compute Jacobian elements - // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption - // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. - - auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; - - result.J00 = 1.0; - result.J0g.fillin(cscale); - result.Jg0 = 1.0 / c_v * dEg_dT; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] <= 0.0) { - result.Jgg[g] = -std::numeric_limits::infinity(); - } else { - result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; - } - } - - return result; -} - -// Compute the Jacobian of energy update equations for the gas-dust-radiation system. The result is a struct containing the following elements: -// J00: (0, 0) component of the Jacobian matrix. = d F0 / d Egas -// F0: (0) component of the residual. = Egas residual -// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 -// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g -// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d Egas -// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g -// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual -template -AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDust( - double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, - quokka::valarray const &Src, double coeff_n, quokka::valarray const &tau, double c_v, double /*lambda_gd_time_dt*/, - quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult -{ - JacobianResult result; - - const double cscale = c_light_ / c_hat_; - - result.F0 = Egas_diff; - result.Fg = Erad_diff - (Rvec + Src); - result.Fg_abs_sum = 0.0; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] > 0.0) { - result.Fg_abs_sum += std::abs(result.Fg[g]); - result.F0 += cscale * Rvec[g]; - } - } - - // const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivativeMultiGroup(T_d, radBoundaries_g_copy); - AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); - - // compute Jacobian elements - // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption - // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. - - auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; - - result.J00 = 1.0; - result.J0g.fillin(cscale); - const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); - // const double coeff_n = dt * dustGasCoeff_local * num_den * num_den / cscale; - dEg_dT *= d_Td_d_T; - const double dTd_dRg = -1.0 / (coeff_n * std::sqrt(T_gas)); - const auto rg = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg; - result.Jg0 = 1.0 / c_v * dEg_dT - 1.0 / cscale * rg * result.J00; - // Note that Fg is modified here, but it does not change Fg_abs_sum, which is used to check the convergence. - result.Fg = result.Fg - 1.0 / cscale * rg * result.F0; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] <= 0.0) { - result.Jgg[g] = -std::numeric_limits::infinity(); - } else { - result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; - } - } - - return result; -} - -// Compute the Jacobian of energy update equations for the gas-dust-radiation system with gas and dust decoupled. The result is a struct containing the -// following elements: J00: (0, 0) component of the Jacobian matrix. = d F0 / d T_d F0: (0) component of the residual. = sum_g R_g - lambda_gd_time_dt -// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 -// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g -// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d T_d -// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g -// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual -template -AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDustDecoupled( - double /*T_gas*/, double /*T_d*/, double /*Egas_diff*/, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, - quokka::valarray const &Src, double /*coeff_n*/, quokka::valarray const &tau, double /*c_v*/, double lambda_gd_time_dt, - quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult -{ - JacobianResult result; - - result.F0 = -lambda_gd_time_dt; - result.Fg = Erad_diff - (Rvec + Src); - result.Fg_abs_sum = 0.0; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] > 0.0) { - result.F0 += Rvec[g]; - result.Fg_abs_sum += std::abs(result.Fg[g]); - } - } - - // compute Jacobian elements - // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption - // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. - - auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; - - result.J00 = 0.0; - result.J0g.fillin(1.0); - result.Jg0 = dEg_dT; - for (int g = 0; g < nGroups_; ++g) { - if (tau[g] <= 0.0) { - result.Jgg[g] = -std::numeric_limits::infinity(); - } else { - result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; - } - } - - return result; -} - // Compute kappaE and kappaP based on the opacity model. The result is stored in the last five arguments: alpha_P, alpha_E, kappaP, kappaE, and kappaPoverE. template AMREX_GPU_DEVICE auto RadSystem::ComputeModelDependentKappaEAndKappaP( @@ -250,14 +94,61 @@ RadSystem::ComputeModelDependentKappaFAndDeltaTerms(double const T, d } } +// Compute the Jacobian of energy update equations for the gas-radiation system. The result is a struct containing the following elements: +// J00: (0, 0) component of the Jacobian matrix. = d F0 / d Egas +// F0: (0) component of the residual. = Egas residual +// Fg_abs_sum: sum of the absolute values of the each component of Fg that has tau(g) > 0 +// J0g: (0, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d F0 / d R_g +// Jg0: (g, 0) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d Egas +// Jgg: (g, g) components of the Jacobian matrix, g = 1, 2, ..., nGroups. = d Fg / d R_g +// Fg: (g) components of the residual, g = 1, 2, ..., nGroups. = Erad residual template -template -AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( - double const Egas0, quokka::valarray const &Erad0Vec, double const rho, double const T_d0, int const dust_model, double const coeff_n, - double const lambda_gd_times_dt, double const dt, amrex::GpuArray const &massScalars, int const n_outer_iter, - quokka::valarray const &work, quokka::valarray const &vel_times_F, quokka::valarray const &Src, - amrex::GpuArray const &rad_boundaries, JacobianFunc ComputeJacobian, int *p_iteration_counter, int *p_iteration_failure_counter) - -> NewtonIterationResult +AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGas(double /*T_d*/, double Egas_diff, quokka::valarray const &Erad_diff, + quokka::valarray const &Rvec, quokka::valarray const &Src, + quokka::valarray const &tau, double c_v, + quokka::valarray const &kappaPoverE, + quokka::valarray const &d_fourpiboverc_d_t) -> JacobianResult +{ + JacobianResult result; + + const double cscale = c_light_ / c_hat_; + + result.F0 = Egas_diff; + result.Fg = Erad_diff - (Rvec + Src); + result.Fg_abs_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] > 0.0) { + result.Fg_abs_sum += std::abs(result.Fg[g]); + result.F0 += cscale * Rvec[g]; + } + } + + // compute Jacobian elements + // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption + // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. + + auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; + + result.J00 = 1.0; + result.J0g.fillin(cscale); + result.Jg0 = 1.0 / c_v * dEg_dT; + for (int g = 0; g < nGroups_; ++g) { + if (tau[g] <= 0.0) { + result.Jgg[g] = -std::numeric_limits::infinity(); + } else { + result.Jgg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; + } + } + + return result; +} + +template +AMREX_GPU_DEVICE auto RadSystem::SolveGasRadiationEnergyExchange( + double const Egas0, quokka::valarray const &Erad0Vec, double const rho, double const dt, + amrex::GpuArray const &massScalars, int const n_outer_iter, quokka::valarray const &work, + quokka::valarray const &vel_times_F, quokka::valarray const &Src, + amrex::GpuArray const &rad_boundaries, int *p_iteration_counter, int *p_iteration_failure_counter) -> NewtonIterationResult { // 1. Compute energy exchange @@ -283,16 +174,10 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( const double cscale = c / chat; // const double Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); - double Etot0 = NAN; - if (dust_model == 0 || dust_model == 1) { - Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); - } else { - // for dust_model == 2 (decoupled gas and dust), Egas0 is not involved in the iteration - Etot0 = std::abs(lambda_gd_times_dt) + (sum(Erad0Vec) + sum(Src)); - } + double Etot0 = Egas0 + cscale * (sum(Erad0Vec) + sum(Src)); double T_gas = NAN; - double T_d = NAN; + double T_d = NAN; // a dummy dust temperature, T_d = T_gas for gas-only model double delta_x = NAN; quokka::valarray delta_R{}; quokka::valarray F_D{}; @@ -305,7 +190,7 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( amrex::GpuArray, 2> kappa_expo_and_lower_value{}; OpacityTerms opacity_terms{}; - // fill kappa_expo_and_lower_value with NAN to avoid mistakes caused by uninitialized values + // fill kappa_expo_and_lower_value with NAN to get warned when there are uninitialized values for (int i = 0; i < 2; ++i) { for (int j = 0; j < nGroups_ + 1; ++j) { kappa_expo_and_lower_value[i][j] = NAN; @@ -338,10 +223,6 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( double Egas_guess = Egas0; auto EradVec_guess = Erad0Vec; - if (dust_model == 2) { - Egas_guess = Egas0 - cscale * lambda_gd_times_dt; // update Egas_guess once for all - } - const double resid_tol = 1.0e-11; // 1.0e-15; const int maxIter = 100; int n = 0; @@ -351,24 +232,7 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); AMREX_ASSERT(T_gas >= 0.); - - if (dust_model == 0) { - T_d = T_gas; - } else if (dust_model == 1) { - if (n == 0) { - T_d = T_d0; - } else { - T_d = T_gas - sum(Rvec) / (coeff_n * std::sqrt(T_gas)); - } - } else if (dust_model == 2) { - if (n == 0) { - T_d = T_d0; - } - } - AMREX_ASSERT_WITH_MESSAGE(T_d >= 0., "Dust temperature is negative!"); - if (T_d < 0.0) { - amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT - } + T_d = T_gas; // 2. Compute kappaP and kappaE at dust temperature @@ -379,7 +243,7 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( if (n == 0) { // Compute kappaF and the delta_nu_kappa_B term. kappaF is used to compute the work term. - // Only the last two arguments (kappaFVec, delta_nu_kappa_B_at_edge) are modified in this function. + // Will update opacity_terms in place ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place } @@ -441,10 +305,8 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( const auto Egas_diff = Egas_guess - Egas0; const auto Erad_diff = EradVec_guess - Erad0Vec; - JacobianResult jacobian; - jacobian = ComputeJacobian(T_gas, T_d, Egas_diff, Erad_diff, Rvec, Src, coeff_n, tau, c_v, lambda_gd_times_dt, opacity_terms.kappaPoverE, - d_fourpiboverc_d_t); + auto jacobian = ComputeJacobianForGas(T_d, Egas_diff, Erad_diff, Rvec, Src, tau, c_v, opacity_terms.kappaPoverE, d_fourpiboverc_d_t); if constexpr (use_D_as_base) { jacobian.J0g = jacobian.J0g * tau0; @@ -461,26 +323,22 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( if (n == 0) { std::cout << "Egas0 = " << Egas0 << ", Erad0Vec = ["; for (int g = 0; g < nGroups_; ++g) { - std::cout << Erad0Vec[g]; - if (g < nGroups_ - 1) { std::cout << ", "; } + std::cout << Erad0Vec[g] << ", "; } std::cout << "], tau0 = ["; for (int g = 0; g < nGroups_; ++g) { - std::cout << tau0[g]; - if (g < nGroups_ - 1) { std::cout << ", "; } + std::cout << tau0[g] << ", "; } std::cout << "]"; std::cout << "; C_V = " << c_v << ", a_rad = " << radiation_constant_ << ", coeff_n = " << coeff_n << "\n"; } else if (n >= 0) { std::cout << "n = " << n << ", Egas_guess = " << Egas_guess << ", EradVec_guess = ["; for (int g = 0; g < nGroups_; ++g) { - std::cout << EradVec_guess[g]; - if (g < nGroups_ - 1) { std::cout << ", "; } + std::cout << EradVec_guess[g] << ", "; } std::cout << "], tau = ["; for (int g = 0; g < nGroups_; ++g) { - std::cout << tau[g]; - if (g < nGroups_ - 1) { std::cout << ", "; } + std::cout << tau[g] << ", "; } std::cout << "]"; std::cout << ", F_G = " << jacobian.F0 << ", F_D_abs_sum = " << jacobian.Fg_abs_sum << ", Etot0 = " << Etot0 << "\n"; @@ -495,21 +353,16 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( // Update independent variables (Egas_guess, Rvec) // enable_dE_constrain is used to prevent the gas temperature from dropping/increasing below/above the radiation // temperature - if (dust_model == 2) { - T_d += delta_x; - Rvec += delta_R; + const double T_rad = std::sqrt(std::sqrt(sum(EradVec_guess) / radiation_constant_)); + if (enable_dE_constrain && delta_x / c_v > std::max(T_gas, T_rad)) { + Egas_guess = quokka::EOS::ComputeEintFromTgas(rho, T_rad); + // Rvec.fillin(0.0); } else { - const double T_rad = std::sqrt(std::sqrt(sum(EradVec_guess) / radiation_constant_)); - if (enable_dE_constrain && delta_x / c_v > std::max(T_gas, T_rad)) { - Egas_guess = quokka::EOS::ComputeEintFromTgas(rho, T_rad); - // Rvec.fillin(0.0); + Egas_guess += delta_x; + if constexpr (use_D_as_base) { + Rvec += tau0 * delta_R; } else { - Egas_guess += delta_x; - if constexpr (use_D_as_base) { - Rvec += tau0 * delta_R; - } else { - Rvec += delta_R; - } + Rvec += delta_R; } } @@ -535,7 +388,7 @@ AMREX_GPU_DEVICE auto RadSystem::SolveMatterRadiationEnergyExchange( if (n > 0) { // calculate kappaF since the temperature has changed - // Only the last two arguments (kappaFVec, delta_nu_kappa_B_at_edge) are modified in this function. + // Will update opacity_terms in place ComputeModelDependentKappaFAndDeltaTerms(T_d, rho, rad_boundaries, fourPiBoverC, opacity_terms); // update opacity_terms in place } @@ -827,29 +680,6 @@ void RadSystem::AddSourceTermsMultiGroup(array_t &consVar, arrayconst if constexpr (gamma_ != 1.0) { - // 1.1. If enable_dust_model, determine if or not to use the weak-coupling approximation - - int dust_model = 0; - double T_d0 = NAN; - double lambda_gd_times_dt = NAN; - if constexpr (enable_dust_gas_thermal_coupling_model_) { - const double T_gas0 = quokka::EOS::ComputeTgasFromEint(rho, Egas0, massScalars); - AMREX_ASSERT(T_gas0 >= 0.); - T_d0 = ComputeDustTemperatureBateKeto(T_gas0, T_gas0, rho, Erad0Vec, coeff_n, dt, NAN, 0, radBoundaries_g_copy); - AMREX_ASSERT_WITH_MESSAGE(T_d0 >= 0., "Dust temperature is negative!"); - if (T_d0 < 0.0) { - amrex::Gpu::Atomic::Add(&p_iteration_failure_counter[1], 1); // NOLINT - } - - const double max_Gamma_gd = coeff_n * std::max(std::sqrt(T_gas0) * T_gas0, std::sqrt(T_d0) * T_d0); - if (cscale * max_Gamma_gd < ISM_Traits::gas_dust_coupling_threshold * Egas0) { - dust_model = 2; - lambda_gd_times_dt = coeff_n * std::sqrt(T_gas0) * (T_gas0 - T_d0); - } else { - dust_model = 1; - } - } - // 1.2. Compute a term required to calculate the work. This is only required in the first outer loop. quokka::valarray vel_times_F{}; @@ -869,19 +699,25 @@ void RadSystem::AddSourceTermsMultiGroup(array_t &consVar, arrayconst // the work term. if constexpr (!enable_dust_gas_thermal_coupling_model_) { - updated_energy = SolveMatterRadiationEnergyExchange( - Egas0, Erad0Vec, rho, T_d0, dust_model, coeff_n, lambda_gd_times_dt, dt, massScalars, iter, work, vel_times_F, Src, - radBoundaries_g_copy, &ComputeJacobianForGas, p_iteration_counter_local, p_iteration_failure_counter_local); + // gas + radiation + updated_energy = + SolveGasRadiationEnergyExchange(Egas0, Erad0Vec, rho, dt, massScalars, iter, work, vel_times_F, Src, + radBoundaries_g_copy, p_iteration_counter_local, p_iteration_failure_counter_local); } else { - auto ComputeJacobian = (dust_model == 1) ? &ComputeJacobianForGasAndDust : &ComputeJacobianForGasAndDustDecoupled; - - updated_energy = SolveMatterRadiationEnergyExchange( - Egas0, Erad0Vec, rho, T_d0, dust_model, coeff_n, lambda_gd_times_dt, dt, massScalars, iter, work, vel_times_F, Src, - radBoundaries_g_copy, ComputeJacobian, p_iteration_counter_local, p_iteration_failure_counter_local); + if constexpr (!enable_photoelectric_heating_) { + // gas + radiation + dust + updated_energy = SolveGasDustRadiationEnergyExchange( + Egas0, Erad0Vec, rho, coeff_n, dt, massScalars, iter, work, vel_times_F, Src, radBoundaries_g_copy, + p_iteration_counter_local, p_iteration_failure_counter_local); + } else { + // gas + radiation + dust + photoelectric heating + updated_energy = SolveGasDustRadiationEnergyExchangeWithPE( + Egas0, Erad0Vec, rho, coeff_n, dt, massScalars, iter, work, vel_times_F, Src, radBoundaries_g_copy, + p_iteration_counter_local, p_iteration_failure_counter_local); + } } Egas_guess = updated_energy.Egas; - EradVec_guess = updated_energy.EradVec; // copy work to work_prev for (int g = 0; g < nGroups_; ++g) { @@ -964,7 +800,6 @@ void RadSystem::AddSourceTermsMultiGroup(array_t &consVar, arrayconst consNew(i, j, k, gasInternalEnergy_index) = Egas_guess; consNew(i, j, k, gasEnergy_index) = ComputeEgasFromEint(rho, x1GasMom1, x2GasMom1, x3GasMom1, Egas_guess); } else { - amrex::ignore_unused(EradVec_guess); amrex::ignore_unused(Egas_guess); amrex::ignore_unused(Egas0); amrex::ignore_unused(Etot0); diff --git a/src/util/valarray.hpp b/src/util/valarray.hpp index ddfeb44c1..9b01a3a84 100644 --- a/src/util/valarray.hpp +++ b/src/util/valarray.hpp @@ -86,6 +86,26 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator+(quokka::valarray c return sum; } +// array + scalar +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator+(quokka::valarray const &v, T const &scalar) -> quokka::valarray +{ + quokka::valarray scalarsum; + for (size_t i = 0; i < v.size(); ++i) { + scalarsum[i] = v[i] + scalar; + } + return scalarsum; +} + +// scalar + array +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator+(T const &scalar, quokka::valarray const &v) -> quokka::valarray +{ + quokka::valarray scalarsum; + for (size_t i = 0; i < v.size(); ++i) { + scalarsum[i] = scalar + v[i]; + } + return scalarsum; +} + // array - array template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator-(quokka::valarray const &a, quokka::valarray const &b) -> quokka::valarray @@ -155,16 +175,6 @@ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void opera } } -// array + scalar -template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator+(quokka::valarray const &v, T const &scalar) -> quokka::valarray -{ - quokka::valarray scalarsum; - for (size_t i = 0; i < v.size(); ++i) { - scalarsum[i] = v[i] + scalar; - } - return scalarsum; -} - // array / scalar template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator/(quokka::valarray const &v, T const &scalar) -> quokka::valarray { diff --git a/tests/RadMarshakDust.in b/tests/RadMarshakDust.in index 69587d670..6a0dc52ee 100644 --- a/tests/RadMarshakDust.in +++ b/tests/RadMarshakDust.in @@ -21,7 +21,7 @@ amr.blocking_factor = 4 # grid size must be divisible by this # Radiation # ***************************************************************** radiation.print_iteration_counts = 1 -radiation.dust_gas_interaction_coeff = 1e-10 +radiation.dust_gas_interaction_coeff = 1e-2 # ***************************************************************** # Problem diff --git a/tests/RadMarshakDustPEcoupled.in b/tests/RadMarshakDustPEcoupled.in new file mode 100644 index 000000000..d7b394594 --- /dev/null +++ b/tests/RadMarshakDustPEcoupled.in @@ -0,0 +1,37 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 0 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 256 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this + +# ***************************************************************** +# Radiation +# ***************************************************************** +radiation.print_iteration_counts = 1 +radiation.dust_gas_interaction_coeff = 1e20 + +# ***************************************************************** +# Problem +# ***************************************************************** +problem.kappa1 = 1e-20 +problem.kappa2 = 1e-20 + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 0 +stop_time = 0.5 +#checkpoint_interval = 1 +#restartfile = chk00006 diff --git a/tests/RadMarshakDustPEdecoupled.in b/tests/RadMarshakDustPEdecoupled.in new file mode 100644 index 000000000..9b7f95cb6 --- /dev/null +++ b/tests/RadMarshakDustPEdecoupled.in @@ -0,0 +1,37 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 0 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 256 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this + +# ***************************************************************** +# Radiation +# ***************************************************************** +radiation.print_iteration_counts = 1 +radiation.dust_gas_interaction_coeff = 1e-20 + +# ***************************************************************** +# Problem +# ***************************************************************** +problem.kappa1 = 1e-20 +problem.kappa2 = 1e-20 + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 0 +stop_time = 0.5 +#checkpoint_interval = 1 +#restartfile = chk00006