From 76eff2f67c70e6804ce86f626d56b24a648e7884 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Sat, 31 Aug 2024 04:13:19 +1000 Subject: [PATCH] improve convergence rate of the radiation Newton-Raphson solver (#720) ### Description Changes: 1. Added a trick to make the Newton-Raphson iteration converge faster in some situations: when $\Delta T > \max(T_{\rm gas}, T_{\rm rad})$, set $T_{\rm gas}$ to $T_{\rm rad}$ and $R$ to 0. The motivation for doing this trick is explained below. 2. Added `radiation.print_iteration_counts` parameter to enable counting the mean and max number of Newton-Raphson iterations in each radiation step. When set to `true`, the log file would look like: ```log Coarse STEP 1 at t = 0 (0%) starts ... time_subcycle = 0, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1 Coarse STEP 2 at t = 0.0046875 (0.9375%) starts ... time_subcycle = 0.0046875, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1 ``` 3. Defined a new radiation problem, `RadMarshakDust`. Initially, the gas is set with a density and temperature of 1, and the initial radiation energy density is 0. Radiation with a temperature $T_r = 0.01$ streams in from the left boundary. Before the dust model was introduced, we assumed gas and dust were perfectly coupled, and the gas would be heated by radiation where it could penetrate, and in regions without radiation, the gas temperature would cool to 0.003159 (solution of $1 = a_R T^4 + T$). In this test, however, we turned on the dust model and set the dust-gas interaction coefficient to a very small number. As a result, the dust and gas are decoupled: the dust temperature will reach the radiation temperature and stay in equilibrium, while the gas remains at its initial temperature. In a future PR, I will extend this test by incorporating a multigroup model. This will demonstrate how FUV radiation from the left boundary is absorbed by the dust, which then re-radiates in IR. ### Related issues None. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [ ] I have added a link to any related issues see (see above). - [x] I have read the [[Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md)](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [x] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [x] I have requested a reviewer for this PR. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- src/QuokkaSimulation.hpp | 46 +++- src/problems/CMakeLists.txt | 1 + src/problems/RadMarshakDust/CMakeLists.txt | 7 + .../test_radiation_marshak_dust.cpp | 243 ++++++++++++++++++ .../test_radiation_marshak_dust.hpp | 24 ++ src/radiation/radiation_system.hpp | 68 +++-- tests/RadMarshakDust.in | 29 +++ 7 files changed, 386 insertions(+), 32 deletions(-) create mode 100644 src/problems/RadMarshakDust/CMakeLists.txt create mode 100644 src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp create mode 100644 src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp create mode 100644 tests/RadMarshakDust.in 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