diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 9c72f170b..2985e8a54 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -131,6 +131,7 @@ template class QuokkaSimulation : public AMRSimulation class QuokkaSimulation : public AMRSimulation const &stateNew, const amrex::Box &indexRange, amrex::Real time, double dt, int stage, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, - amrex::GpuArray const &prob_hi, int *num_failed_coupling, int *num_failed_dust, - int *p_num_failed_outer_ite); + amrex::GpuArray const &prob_hi, int *p_iteration_counter, int *num_failed_coupling, + int *num_failed_dust, int *p_num_failed_outer_ite); auto computeRadiationFluxes(amrex::Array4 const &consVar, const amrex::Box &indexRange, int nvars, amrex::GpuArray dx) @@ -390,6 +391,7 @@ template void QuokkaSimulation::readParmParse() rpp.query("reconstruction_order", radiationReconstructionOrder_); rpp.query("cfl", radiationCflNumber_); rpp.query("dust_gas_interaction_coeff", dustGasInteractionCoeff_); + rpp.query("print_iteration_counts", print_rad_counter_); } } @@ -1619,9 +1621,11 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real amrex::Gpu::Buffer num_failed_coupling({0}); amrex::Gpu::Buffer num_failed_dust({0}); amrex::Gpu::Buffer num_failed_outer({0}); + amrex::Gpu::Buffer iteration_counter({0, 0, 0}); int *p_num_failed_coupling = num_failed_coupling.data(); int *p_num_failed_dust = num_failed_dust.data(); int *p_num_failed_outer = num_failed_outer.data(); + int *p_iteration_counter = iteration_counter.data(); if constexpr (IMEX_a22 > 0.0) { // matter-radiation exchange source terms of stage 1 @@ -1634,8 +1638,8 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real // update state_new_cc_[lev] in place (updates both radiation and hydro vars) // Note that only a fraction (IMEX_a32) of the matter-radiation exchange source terms are added to hydro. This ensures that the // hydro properties get to t + IMEX_a32 dt in terms of matter-radiation exchange. - operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi, p_num_failed_coupling, - p_num_failed_dust, p_num_failed_outer); + operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi, p_iteration_counter, + p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer); } } @@ -1652,8 +1656,30 @@ void QuokkaSimulation::subcycleRadiationAtLevel(int lev, amrex::Real auto const &prob_lo = geom[lev].ProbLoArray(); auto const &prob_hi = geom[lev].ProbHiArray(); // update state_new_cc_[lev] in place (updates both radiation and hydro vars) - operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi, p_num_failed_coupling, - p_num_failed_dust, p_num_failed_outer); + operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi, p_iteration_counter, + p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer); + } + + 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 + + amrex::ParallelDescriptor::ReduceLongSum(global_solver_count); + amrex::ParallelDescriptor::ReduceLongSum(global_iteration_sum); + amrex::ParallelDescriptor::ReduceIntMax(global_iteration_max); + + const auto n_cells = CountCells(lev); + const double global_iteration_mean = 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 + + if (amrex::ParallelDescriptor::IOProcessor()) { + amrex::Print() << "time_subcycle = " << time_subcycle << ", total number of cells updated is " << CountCells(lev) + << ", 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"; + } } const int nf_coupling = *(num_failed_coupling.copyToHost()); @@ -1820,8 +1846,8 @@ template void QuokkaSimulation::operatorSplitSourceTerms(amrex::Array4 const &stateNew, const amrex::Box &indexRange, const amrex::Real time, const double dt, const int stage, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, - amrex::GpuArray const &prob_hi, int *p_num_failed_coupling, - int *p_num_failed_dust, int *p_num_failed_outer_ite) + amrex::GpuArray const &prob_hi, int *p_iteration_counter, + int *p_num_failed_coupling, int *p_num_failed_dust, int *p_num_failed_outer_ite) { amrex::FArrayBox radEnergySource(indexRange, Physics_Traits::nGroups, amrex::The_Async_Arena()); // cell-centered scalar @@ -1832,8 +1858,8 @@ void QuokkaSimulation::operatorSplitSourceTerms(amrex::Array4::SetRadEnergySource(radEnergySource.array(), indexRange, dx, prob_lo, prob_hi, time + dt); // cell-centered source terms - RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_, p_num_failed_coupling, - p_num_failed_dust, p_num_failed_outer_ite); + RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_, p_iteration_counter, + p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer_ite); } template diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index 24c827be7..2c6e32528 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -35,6 +35,7 @@ add_subdirectory(RadTophat) add_subdirectory(RadTube) add_subdirectory(RadDust) add_subdirectory(RadDustMG) +add_subdirectory(RadMarshakDust) add_subdirectory(RadhydroShell) add_subdirectory(RadhydroShock) diff --git a/src/problems/RadMarshakDust/CMakeLists.txt b/src/problems/RadMarshakDust/CMakeLists.txt new file mode 100644 index 000000000..be399838f --- /dev/null +++ b/src/problems/RadMarshakDust/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_radiation_marshak_dust test_radiation_marshak_dust.cpp ../../util/fextract.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radiation_marshak_dust) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadiationMarshakDust COMMAND test_radiation_marshak_dust RadMarshakDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp new file mode 100644 index 000000000..5213f2e24 --- /dev/null +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp @@ -0,0 +1,243 @@ +//============================================================================== +// 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.cpp +/// \brief Defines a test Marshak wave problem with weak coupling between dust and gas. +/// + +#include "test_radiation_marshak_dust.hpp" +#include "AMReX.H" +#include "QuokkaSimulation.hpp" +#include "util/fextract.hpp" +#include "util/valarray.hpp" + +struct StreamingProblem { +}; + +constexpr double c = 1.0; // speed of light +constexpr double chat = 1.0; // reduced speed of light +constexpr double kappa0 = 10.0; // opacity +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.0e10; +constexpr double erad_floor = 1.0e-10; +constexpr double initial_Erad = erad_floor; +constexpr double T_rad_L = 1.0e-2; // so EradL = 1e2 +constexpr double EradL = a_rad * T_rad_L * T_rad_L * T_rad_L * T_rad_L; +// constexpr double T_end_exact = 0.0031597766719577; // dust off; solution of 1 == a_rad * T^4 + T +constexpr double T_end_exact = initial_T; // dust on + +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 = 1; // number of radiation groups +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr int beta_order = 0; + static constexpr bool enable_dust_gas_thermal_coupling_model = true; +}; + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real +{ + return kappa0; +} + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real +{ + return kappa0; +} + +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 Erad0 = initial_Erad; + 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) = Erad0; + 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(); + + if (i < lo[0]) { + // streaming inflow boundary + const double Erad = EradL; + const double Frad = c * Erad; + + // 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) = Erad; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frad; + 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; + + // 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 erad(nx); + std::vector T(nx); + std::vector T_exact(nx); + std::vector xs(nx); + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + double erad_sim = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + erad_sim += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + erad.at(i) = erad_sim; + 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; + } + + double err_norm = 0.; + double sol_norm = 0.; + for (int i = 0; i < nx; ++i) { + err_norm += std::abs(T[i] - T_exact[i]); + sol_norm += std::abs(T_exact[i]); + } + + const double rel_err_norm = err_norm / sol_norm; + const double rel_err_tol = 0.02; + int status = 1; + if (rel_err_norm < rel_err_tol) { + status = 0; + } + amrex::Print() << "Relative L1 norm = " << rel_err_norm << std::endl; + +#ifdef HAVE_PYTHON + // Plot results + matplotlibcpp::clf(); + std::map plot_args; + plot_args["label"] = "numerical solution"; + matplotlibcpp::plot(xs, erad, plot_args); + matplotlibcpp::xlabel("x"); + matplotlibcpp::ylabel("E_rad"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("t = {:f}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radiation_marshak_dust_Erad.pdf"); + + // plot temperature + matplotlibcpp::clf(); + matplotlibcpp::ylim(0.0, 1.1); + matplotlibcpp::plot(xs, T, plot_args); + matplotlibcpp::xlabel("x"); + matplotlibcpp::ylabel("Temperature"); + matplotlibcpp::title(fmt::format("t = {:f}", sim.tNew_[0])); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radiation_marshak_dust_temperature.pdf"); +#endif // HAVE_PYTHON + + // Cleanup and exit + amrex::Print() << "Finished." << std::endl; + return status; +} diff --git a/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp b/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp new file mode 100644 index 000000000..7ece497a8 --- /dev/null +++ b/src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp @@ -0,0 +1,24 @@ +#ifndef TEST_RADIATION_MARSHAK_DUST_HPP_ // NOLINT +#define TEST_RADIATION_MARSHAK_DUST_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_system.hpp" + +// function definitions + +#endif // TEST_RADIATION_MARSHAK_DUST_HPP_ diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 99fe99393..f9aaa770d 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -33,18 +33,16 @@ using Real = amrex::Real; -// Hyper parameters of the radiation solver - +// Hyper parameters for the radiation solver static constexpr bool include_delta_B = true; static constexpr bool use_diffuse_flux_mean_opacity = true; static constexpr bool special_edge_bin_slopes = false; // Use 2 and -4 as the slopes for the first and last bins, respectively static constexpr bool force_rad_floor_in_iteration = false; // force radiation energy density to be positive (and above the floor value) in the Newton iteration +static constexpr bool include_work_term_in_source = true; -static constexpr bool use_diffuse_flux_mean_opacity_incl_kappaE = true; static const int max_ite_to_update_alpha_E = 5; // Apply to the PPL_opacity_full_spectrum only. Only update alpha_E for the first max_ite_to_update_alpha_E // iterations of the Newton iteration - -static constexpr bool include_work_term_in_source = true; +static constexpr bool enable_dE_constrain = true; static constexpr bool use_D_as_base = false; static const bool PPL_free_slope_st_total = false; // PPL with free slopes for all, but subject to the constraint sum_g alpha_g B_g = - sum_g B_g. Not working // well -- Newton iteration convergence issue. @@ -205,7 +203,7 @@ template class RadSystem : public HyperbolicSystem::ComputeDustTemperature(double c template void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt_radiation, - const int stage, double dustGasCoeff, int *p_num_failed_coupling, int *p_num_failed_dust, int *p_num_failed_outer_ite) + const int stage, double dustGasCoeff, int *p_iteration_counter, int *p_num_failed_coupling, int *p_num_failed_dust, + int *p_num_failed_outer_ite) { arrayconst_t &consPrev = consVar; // make read-only array_t &consNew = consVar; @@ -1286,6 +1285,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne auto p_num_failed_coupling_local = p_num_failed_coupling; auto p_num_failed_dust_local = p_num_failed_dust; auto p_num_failed_outer_local = p_num_failed_outer_ite; + auto p_iteration_counter_local = p_iteration_counter; const double c = c_light_; const double chat = c_hat_; @@ -1447,7 +1447,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray F_D{}; const double resid_tol = 1.0e-11; // 1.0e-15; - const int maxIter = 50; + const int maxIter = enable_dust_gas_thermal_coupling_model_ ? 100 : 50; int n = 0; for (; n < maxIter; ++n) { T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); @@ -1609,13 +1609,6 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } - if (n > 200) { -#ifndef NDEBUG - std::cout << "n = " << n << ", F_G = " << F_G << ", F_D_abs_sum = " << F_D_abs_sum - << ", F_D_abs_sum / Etot0 = " << F_D_abs_sum / Etot0 << std::endl; -#endif - } - // check relative convergence of the residuals if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * F_D_abs_sum / Etot0 < resid_tol)) { break; @@ -1623,6 +1616,18 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T +#if 0 + // For debugging: print (Egas0, Erad0Vec, tau0), which defines the initial condition for a Newton-Raphson iteration + if (n == maxIter - 10) { + std::cout << "Egas0 = " << Egas0 << ", Erad0Vec = " << Erad0Vec[0] << ", tau0 = " << tau0[0] + << "; C_V = " << c_v << ", a_rad = " << radiation_constant_ << std::endl; + } else if (n >= maxIter - 10) { + std::cout << "n = " << n << ", Egas_guess = " << Egas_guess << ", EradVec_guess = " << EradVec_guess[0] + << ", tau = " << tau[0]; + std::cout << ", F_G = " << F_G << ", F_D_abs_sum = " << F_D_abs_sum << ", Etot0 = " << Etot0 << std::endl; + } +#endif + const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivative(T_d, radBoundaries_g_copy); AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); @@ -1675,11 +1680,26 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(!std::isnan(deltaEgas)); AMREX_ASSERT(!deltaD.hasnan()); - Egas_guess += deltaEgas; - if constexpr (use_D_as_base) { - Rvec += tau0 * deltaD; + if (!enable_dE_constrain) { + Egas_guess += deltaEgas; + if constexpr (use_D_as_base) { + Rvec += tau0 * deltaD; + } else { + Rvec += deltaD; + } } else { - Rvec += deltaD; + const double T_rad = std::pow(sum(EradVec_guess) / radiation_constant_, 0.25); + if (deltaEgas / c_v > std::max(T_gas, T_rad)) { + Egas_guess = quokka::EOS::ComputeEintFromTgas(rho, T_rad); + Rvec.fillin(0.0); + } else { + Egas_guess += deltaEgas; + if constexpr (use_D_as_base) { + Rvec += tau0 * deltaD; + } else { + Rvec += deltaD; + } + } } // check relative and absolute convergence of E_r @@ -1693,6 +1713,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::Gpu::Atomic::Add(p_num_failed_coupling_local, 1); } + // update iteration counter: (+1, +ite, max(self, ite)) + amrex::Gpu::Atomic::Add(&p_iteration_counter_local[0], 1); // total number of radiation updates + amrex::Gpu::Atomic::Add(&p_iteration_counter_local[1], n + 1); // total number of Newton-Raphson iterations + amrex::Gpu::Atomic::Max(&p_iteration_counter_local[2], n + 1); // maximum number of Newton-Raphson iterations + // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ASSERT(Egas_guess > 0.0); AMREX_ASSERT(min(EradVec_guess) >= 0.0); @@ -1941,9 +1966,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop const double lag_tol = 1.0e-13; - if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) / Etot0 < lag_tol) || - (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || - (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { + if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) < lag_tol * Etot0) || + (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || (sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { break; } } diff --git a/tests/RadMarshakDust.in b/tests/RadMarshakDust.in new file mode 100644 index 000000000..4a8b92648 --- /dev/null +++ b/tests/RadMarshakDust.in @@ -0,0 +1,29 @@ +# ***************************************************************** +# 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 = 64 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 + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 0 +stop_time = 0.5 \ No newline at end of file