diff --git a/docs/references.bib b/docs/references.bib index 3a6bcc5d6..0cc7df648 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -263,3 +263,20 @@ @ARTICLE{Jin_1996 adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@article{Krumholz_2007, + title = {Equations and {{Algorithms}} for {{Mixed-frame Flux-limited Diffusion Radiation Hydrodynamics}}}, + author = {{Krumholz}, Mark R. and {Klein}, Richard I. and {McKee}, Christopher F. and {Bolstad}, John}, + year = {2007}, + month = sep, + journal = {The Astrophysical Journal}, + volume = {667}, + pages = {626--643}, + issn = {0004-637X}, + doi = {10.1086/520791}, + url = {https://ui.adsabs.harvard.edu/abs/2007ApJ...667..626K}, + urldate = {2023-10-18}, + abstract = {We analyze the mixed-frame equations of radiation hydrodynamics under the approximations of flux-limited diffusion and a thermal radiation field and derive the minimal set of evolution equations that includes all terms that are of leading order in any regime of nonrelativistic radiation hydrodynamics. Our equations are accurate to first order in v/c in the static diffusion regime. In contrast, we show that previous lower order derivations of these equations omit leading terms in at least some regimes. In comparison to comoving-frame formulations of radiation hydrodynamics, our equations have the advantage that they manifestly conserve total energy, making them very well suited to numerical simulations, particularly with adaptive meshes. For systems in the static diffusion regime, our analysis also suggests an algorithm that is both simpler and faster than earlier comoving-frame methods. We implement this algorithm in the Orion adaptive mesh refinement code and show that it performs well in a range of test problems.}, + keywords = {Astrophysics,Hydrodynamics,Methods: Numerical,Radiative Transfer}, + annotation = {ADS Bibcode: 2007ApJ...667..626K}, + file = {/Users/cche/Googledrive2/Academic/Zotero/Code/QUOKKA/Krumholz_plus_2007_Equations_and_Algorithms_for_Mixed-frame2.pdf} +} diff --git a/docs/tests/hydro_shuosher.png b/docs/tests/attach/hydro_shuosher.png similarity index 100% rename from docs/tests/hydro_shuosher.png rename to docs/tests/attach/hydro_shuosher.png diff --git a/docs/tests/hydro_sms.png b/docs/tests/attach/hydro_sms.png similarity index 100% rename from docs/tests/hydro_sms.png rename to docs/tests/attach/hydro_sms.png diff --git a/docs/tests/radcoupling.png b/docs/tests/attach/radcoupling.png similarity index 100% rename from docs/tests/radcoupling.png rename to docs/tests/attach/radcoupling.png diff --git a/docs/tests/radcoupling_rsla.png b/docs/tests/attach/radcoupling_rsla.png similarity index 100% rename from docs/tests/radcoupling_rsla.png rename to docs/tests/attach/radcoupling_rsla.png diff --git a/docs/tests/attach/radhydro_pulse_density-1.png b/docs/tests/attach/radhydro_pulse_density-1.png new file mode 100644 index 000000000..73b9c3b41 Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_density-1.png differ diff --git a/docs/tests/attach/radhydro_pulse_density.png b/docs/tests/attach/radhydro_pulse_density.png new file mode 100644 index 000000000..163980b89 Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_density.png differ diff --git a/docs/tests/attach/radhydro_pulse_temperature-1.png b/docs/tests/attach/radhydro_pulse_temperature-1.png new file mode 100644 index 000000000..94a8b8807 Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_temperature-1.png differ diff --git a/docs/tests/attach/radhydro_pulse_temperature.png b/docs/tests/attach/radhydro_pulse_temperature.png new file mode 100644 index 000000000..09167748f Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_temperature.png differ diff --git a/docs/tests/attach/radhydro_pulse_velocity-1.png b/docs/tests/attach/radhydro_pulse_velocity-1.png new file mode 100644 index 000000000..27fa572a7 Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_velocity-1.png differ diff --git a/docs/tests/attach/radhydro_pulse_velocity.png b/docs/tests/attach/radhydro_pulse_velocity.png new file mode 100644 index 000000000..1417764fb Binary files /dev/null and b/docs/tests/attach/radhydro_pulse_velocity.png differ diff --git a/docs/tests/attach/radhydro_uniform_advecting_temperature-nobeta.png b/docs/tests/attach/radhydro_uniform_advecting_temperature-nobeta.png new file mode 100644 index 000000000..50da3eaf4 Binary files /dev/null and b/docs/tests/attach/radhydro_uniform_advecting_temperature-nobeta.png differ diff --git a/docs/tests/attach/radhydro_uniform_advecting_temperature.png b/docs/tests/attach/radhydro_uniform_advecting_temperature.png new file mode 100644 index 000000000..714a6c5d1 Binary files /dev/null and b/docs/tests/attach/radhydro_uniform_advecting_temperature.png differ diff --git a/docs/tests/attach/radhydro_uniform_advecting_velocity-nobeta.png b/docs/tests/attach/radhydro_uniform_advecting_velocity-nobeta.png new file mode 100644 index 000000000..af5e1fb05 Binary files /dev/null and b/docs/tests/attach/radhydro_uniform_advecting_velocity-nobeta.png differ diff --git a/docs/tests/attach/radhydro_uniform_advecting_velocity.png b/docs/tests/attach/radhydro_uniform_advecting_velocity.png new file mode 100644 index 000000000..061c654b8 Binary files /dev/null and b/docs/tests/attach/radhydro_uniform_advecting_velocity.png differ diff --git a/docs/tests/radshock_cgs_temperature.png b/docs/tests/attach/radshock_cgs_temperature.png similarity index 100% rename from docs/tests/radshock_cgs_temperature.png rename to docs/tests/attach/radshock_cgs_temperature.png diff --git a/docs/tests/energy_exchange.rst b/docs/tests/energy_exchange.rst index 655a4c969..b69957e74 100644 --- a/docs/tests/energy_exchange.rst +++ b/docs/tests/energy_exchange.rst @@ -33,13 +33,13 @@ The exact time-dependent solution for the matter temperature :math:`T` is: We show the numerical results below: -.. figure:: radcoupling_rsla.png +.. figure:: attach/radcoupling_rsla.png :alt: A figure showing the radiation temperature and material temperature as a function of time. The radiation temperature and matter temperatures in the reduced speed-of-light approximation, along with the exact solution for the matter temperature. -.. figure:: radcoupling.png +.. figure:: attach/radcoupling.png :alt: A figure showing the radiation temperature and material temperature as a function of time. The radiation temperature and matter temperatures, along with the exact solution for the matter temperature. \ No newline at end of file diff --git a/docs/tests/index.rst b/docs/tests/index.rst index 46b7da098..7b3afa926 100644 --- a/docs/tests/index.rst +++ b/docs/tests/index.rst @@ -10,4 +10,6 @@ Listed here are the test problems that are included with Quokka. *This page is s shu_osher sms energy_exchange + radhydro_uniform_adv + radhydro_pulse \ No newline at end of file diff --git a/docs/tests/radhydro_pulse.rst b/docs/tests/radhydro_pulse.rst new file mode 100644 index 000000000..63f2d8e52 --- /dev/null +++ b/docs/tests/radhydro_pulse.rst @@ -0,0 +1,83 @@ +.. Advecting radiation pulse test + +Advecting radiation pulse test +========================= + +This test demonstrates the code’s ability to deal with the relativistic +correction source terms that arise from the mixed frame formulation of +the RHD moment equations, in a fully-coupled RHD problem. The problems +involve the advection of the a pulse of radiation energy in an optically +thick (:math:`\tau \gg 1`) gas in both static (:math:`\beta \tau \ll 1`) +and dynamic (:math:`\beta \tau \gg 1`) diffusion regimes, with a uniform +background flow velocity :cite:`Krumholz_2007`. + +Parameters +---------- + +Initial condition of the problem in static diffusion regime: + +.. math:: + + + \begin{align} + T = T_0 + (T_1 - T_0) \exp \left( - \frac{x^2}{2 w^2} \right), \\ + w = 24 ~{\rm cm}, T_0 = 10^7 ~{\rm K}, T_1 = 2 \times 10^7 ~{\rm K} \\ + \rho=\rho_0 \frac{T_0}{T}+\frac{a_{\mathrm{R}} \mu}{3 k_{\mathrm{B}}}\left(\frac{T_0^4}{T}-T^3\right) \\ + \rho_0 = 1.2 ~{\rm g~cm^{-3}}, \mu = 2.33 ~m_{\rm H} \\ + \kappa_P=\kappa_R=\kappa = 100 \mathrm{~cm}^2 \mathrm{~g}^{-1} \\ + v = 10 ~{\rm km~s^{-1}} \\ + \tau = \rho \kappa w = 3 \times 10^3, \beta = v/c = 3 \times 10^{-5}, \beta \tau = 9 \times 10^{-2} + \end{align} + +The simulation is run till +:math:`t_{\rm end} = 2 w/v = 4.8 \times 10^{-5} ~{\rm s}`. + +Initial condition of the problem in dynamic diffusion regime: same +parameters as in the static diffusion regime except + +.. math:: + + + \begin{align} + \kappa_P=\kappa_R=\kappa=1000 \mathrm{~cm}^2 \mathrm{~g}^{-1} \\ + v = 1000 ~{\rm km~s^{-1}} \\ + t_{\rm end} = 2 w/v = 1.2 \times 10^{-4} ~{\rm s} \\ + \tau = \rho \kappa w = 3 \times 10^4, \beta = v/c = 3 \times 10^{-3}, \beta \tau = 90 + \end{align} + +Results +------- + +Static diffusion regime: + +.. figure:: attach/radhydro_pulse_temperature-1.png + :alt: radhydro_pulse_temperature-static-diffusion + + radhydro_pulse_temperature-static-diffusion + +.. figure:: attach/radhydro_pulse_density-1.png + :alt: radhydro_pulse_density-static-diffusion + + radhydro_pulse_density-static-diffusion + +.. figure:: attach/radhydro_pulse_velocity-1.png + :alt: radhydro_pulse_velocity-static-diffusion + + radhydro_pulse_velocity-static-diffusion + +Dynamic diffusion regime: + +.. figure:: attach/radhydro_pulse_temperature.png + :alt: radhydro_pulse_temperature-dynamic-diffusion + + radhydro_pulse_temperature-dynamic-diffusion + +.. figure:: attach/radhydro_pulse_density.png + :alt: radhydro_pulse_density-dynamic-diffusion + + radhydro_pulse_density-dynamic-diffusion + +.. figure:: attach/radhydro_pulse_velocity.png + :alt: radhydro_pulse_velocity-dynamic-diffusion + + radhydro_pulse_velocity-dynamic-diffusion diff --git a/docs/tests/radhydro_uniform_adv.rst b/docs/tests/radhydro_uniform_adv.rst new file mode 100644 index 000000000..b52e8278e --- /dev/null +++ b/docs/tests/radhydro_uniform_adv.rst @@ -0,0 +1,95 @@ +Uniform advecting radiation in diffusive limit +============================================== + +In this test, we simulation an advecting uniform gas where radiation and +matter are in thermal equilibrium in the co-moving frame. Following the +Lorentz tranform, the initial radiation energy and flux in the lab frame +to first order in :math:`v/c` are :math:`E_r = a_r T^4` and +:math:`F_r = \frac{4}{3} v E_r`. + +Parameters +---------- + +.. math:: + + + \begin{align} + T_0 = 10^7~{\rm K} \\ + \rho_0 = 1.2 ~{\rm g~cm^{-3}}, \mu = 2.33 ~m_{\rm H} \\ + \kappa_P=\kappa_R=100 \mathrm{~cm}^2 \mathrm{~g}^{-1} \\ + v_{x,0} = 10 ~{\rm km~s^{-1}} \\ + E_{r,0} = a_r T_0^4 \\ + F_{x,0} = \frac{4}{3} v_{x,0} E_{r,0} \\ + t_{\rm end} = 4.8 \times 10^{-5} ~{\rm s} + \end{align} + +Results +------- + +With :math:`O(\beta \tau)` terms: + +.. figure:: attach/radhydro_uniform_advecting_temperature.png + :alt: A figure showing the radiation temperature and material temperature as a function of time. + + The radiation temperature and matter temperatures, along with the exact solution. + +.. figure:: attach/radhydro_uniform_advecting_velocity.png + :alt: A figure showing the radiation velocity and material velocity as a function of time. + + The matter velocity, along with the exact solution. + +Without :math:`O(\beta \tau)` terms: + +.. figure:: attach/radhydro_uniform_advecting_temperature-nobeta.png + :alt: A figure showing the radiation temperature and material temperature as a function of time. + + The radiation temperature and matter temperatures, along with the exact solution. + +.. figure:: attach/radhydro_uniform_advecting_velocity-nobeta.png + :alt: A figure showing the radiation velocity and material velocity as a function of time. + + The matter velocity, along with the exact solution. + + +Physics +------- + +In the transport equation, both the radiation energy and flux are +unchanged because the radiation flux and pressure are uniform. In the +matter-radiation exchange step, the source term is zero since the +radiation and matter are in equilibrium. Finally, the flux is updated +following + +.. math:: + + + \mathbf{F}_{r}^{(t+1)} = \frac{\mathbf{F}_{r}^{(t)} + \Delta t \left[ \rho \kappa_P \left(\frac{4 \pi B}{c}\right) \boldsymbol{v}c + \rho \kappa_F (\boldsymbol{v} :\boldsymbol{P}_r) c \right] }{1+\rho \kappa_{F} {c} \Delta t}. + +With :math:`F_{r}^{(t)} = 4 v E_{r}^{(t)} / 3`, and +:math:`\kappa_P=\kappa_R=\kappa`, we have + +.. math:: + + + \mathbf{F}_{r}^{(t+1)} = \frac{\frac{4}{3} v E_r^{(t)} + \Delta t \left[ \rho \kappa E_r^{(t)} \boldsymbol{v}c + \rho \kappa \boldsymbol{v} (\frac{1}{3}E_r^{(t)}) c \right] }{1+\rho \kappa {c} \Delta t} = \frac{4}{3} v E_r^{(t)} = F_{r}^{(t)} + +Therefore, :math:`F_r` remains constant. This demonstrates that the code +is invariant under Lorentz transformation. + +We can also show that, with the :math:`O(\beta \tau)` terms in the +matter-radiation exchange step, the space-like component of the +radiation four-force vanishes: + +.. math:: + + + \begin{align} + -G &= -\rho \kappa_F \frac{\boldsymbol{F}_r}{c} + \rho \kappa_P\left(\frac{4 \pi B}{c}\right) \frac{\boldsymbol{v}}{c}+\rho \kappa_F \frac{\boldsymbol{v} :\boldsymbol{P}_r}{c} \\ + &= -\rho \kappa \frac{4}{3} E_r v / c + \rho \kappa E_r v / c+ \rho \kappa \frac{1}{3} E_r v / c \\ + &= 0 + \end{align} + +.. |radhydro_uniform_advecting_temperature| image:: attach/radhydro_uniform_advecting_temperature.png +.. |radhydro_uniform_advecting_velocity| image:: attach/radhydro_uniform_advecting_velocity.png +.. |radhydro_uniform_advecting_temperature-nobeta| image:: attach/radhydro_uniform_advecting_temperature-nobeta.png +.. |radhydro_uniform_advecting_velocity-nobeta| image:: attach/radhydro_uniform_advecting_velocity-nobeta.png diff --git a/docs/tests/radshock.rst b/docs/tests/radshock.rst index d651a86a8..b91d561cf 100644 --- a/docs/tests/radshock.rst +++ b/docs/tests/radshock.rst @@ -53,7 +53,7 @@ Since the solution is given assuming radiation diffusion, we set the Eddington f We use the RK2 integrator with a CFL number of 0.2 and a mesh of 256 equally-spaced zones. After 3 shock crossing times, we obtain a solution for the radiation temperature and matter temperature that agrees to better than 0.5% (in relative L1 norm) with the steady-state ODE solution to the radiation hydrodynamics equations: -.. figure:: radshock_cgs_temperature.png +.. figure:: attach/radshock_cgs_temperature.png :alt: A figure showing the radiation temperature and material temperature as a function of position. The radiation temperature is shown in the black solid and dashed lines, with the dashed line showing the semi-analytic solution. The material temperature is shown in the red lines, with the semi-analytic solution shown with the dashed line. \ No newline at end of file diff --git a/docs/tests/shu_osher.rst b/docs/tests/shu_osher.rst index 774d9f7cf..5d0af44c0 100644 --- a/docs/tests/shu_osher.rst +++ b/docs/tests/shu_osher.rst @@ -36,7 +36,7 @@ the interface states in the characteristic variables, as done in §4 of :cite:`S The reference solution is computed using Athena++ with PPM reconstruction in the characteristic variables on a grid of 1600 zones. -.. figure:: hydro_shuosher.png +.. figure:: attach/hydro_shuosher.png :alt: A figure showing the numerical solution to the Shu-Osher test. The density is shown as the solid blue line. There is no exact solution for this problem. \ No newline at end of file diff --git a/docs/tests/sms.rst b/docs/tests/sms.rst index 35606eb99..aa733de7f 100644 --- a/docs/tests/sms.rst +++ b/docs/tests/sms.rst @@ -35,7 +35,7 @@ We use the RK2 integrator with a fixed timestep of :math:`10^{-3}` and a mesh of 100 equally-spaced cells. The contact discontinuity is initially placed at :math:`x=0.5`. -.. figure:: hydro_sms.png +.. figure:: attach/hydro_sms.png :alt: A figure showing the numerical solution to the slow-moving shock test. The density is shown as the solid blue line. The exact solution is the solid orange line. \ No newline at end of file diff --git a/extern/Microphysics b/extern/Microphysics index 657406dca..7de212422 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit 657406dcaf546a5601c8388ef1a2b34a0016edb5 +Subproject commit 7de212422501a88c1c1f536b93217721e3dc6701 diff --git a/extern/amrex b/extern/amrex index 9b733ec45..75571e2dc 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit 9b733ec45cd93a80234a7c98248b6eb4816589d5 +Subproject commit 75571e2dcbf2417529c5ed8e24113580e8e1f3f1 diff --git a/extern/fmt b/extern/fmt index 99b9fbf8e..a5bacf3fe 160000 --- a/extern/fmt +++ b/extern/fmt @@ -1 +1 @@ -Subproject commit 99b9fbf8eff2fc63e639260bd96e22f058f4f818 +Subproject commit a5bacf3fefcbe927b8aafac2ab17d5c6895cabe2 diff --git a/extern/openPMD-api b/extern/openPMD-api index e8e2aeb28..a3fe9b75d 160000 --- a/extern/openPMD-api +++ b/extern/openPMD-api @@ -1 +1 @@ -Subproject commit e8e2aeb28eda917b6b4e21f8d072d55c651b0915 +Subproject commit a3fe9b75d81aa1df569b9fb73ad67b772638961f diff --git a/extern/yaml-cpp b/extern/yaml-cpp index eaf720537..94710bb22 160000 --- a/extern/yaml-cpp +++ b/extern/yaml-cpp @@ -1 +1 @@ -Subproject commit eaf72053724814be3b99d38e292fca5797a57b7b +Subproject commit 94710bb2213cb6793121eaaa4edfb15abcca7af9 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 43fce52a2..fc3a84bf6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -166,6 +166,7 @@ add_subdirectory(RadTube) add_subdirectory(RadhydroShell) add_subdirectory(RadhydroShock) add_subdirectory(RadhydroShockCGS) +add_subdirectory(RadhydroShockMultigroup) add_subdirectory(RadhydroUniformAdvecting) add_subdirectory(RadhydroPulse) add_subdirectory(RadhydroPulseGrey) diff --git a/src/RadhydroPulse/CMakeLists.txt b/src/RadhydroPulse/CMakeLists.txt index 5ea861de9..0cc17b65f 100644 --- a/src/RadhydroPulse/CMakeLists.txt +++ b/src/RadhydroPulse/CMakeLists.txt @@ -1,7 +1,9 @@ -add_executable(test_radhydro_pulse test_radhydro_pulse.cpp ../fextract.cpp ${QuokkaObjSources}) +if (AMReX_SPACEDIM EQUAL 1) + add_executable(test_radhydro_pulse test_radhydro_pulse.cpp ../fextract.cpp ${QuokkaObjSources}) -if(AMReX_GPU_BACKEND MATCHES "CUDA") - setup_target_for_cuda_compilation(test_radhydro_pulse) -endif(AMReX_GPU_BACKEND MATCHES "CUDA") + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_pulse) + endif(AMReX_GPU_BACKEND MATCHES "CUDA") -add_test(NAME RadhydroPulse COMMAND test_radhydro_pulse RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) + add_test(NAME RadhydroPulse COMMAND test_radhydro_pulse RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/RadhydroPulse/test_radhydro_pulse.cpp b/src/RadhydroPulse/test_radhydro_pulse.cpp index 91fecd0dc..5605152cb 100644 --- a/src/RadhydroPulse/test_radhydro_pulse.cpp +++ b/src/RadhydroPulse/test_radhydro_pulse.cpp @@ -29,7 +29,7 @@ constexpr double v0_nonadv = 0.; // non-advecting pulse // static diffusion: tau = 2e3, beta = 3e-5, beta tau = 6e-2 constexpr double kappa0 = 100.; // cm^2 g^-1 constexpr double v0_adv = 1.0e6; // advecting pulse -constexpr double max_time = 4.8e-5; // max_time = 2.0 * width / v1; +constexpr double max_time = 2.4e-5; // max_time = 2.0 * width / v1; // dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 // constexpr double kappa0 = 1000.; // cm^2 g^-1 @@ -350,7 +350,7 @@ auto problem_main() -> int err_norm += std::abs(Tgas2[i] - Trad[i]); sol_norm += std::abs(Trad[i]) * 3.0; } - const double error_tol = 0.001; + const double error_tol = 0.006; const double rel_error = err_norm / sol_norm; amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; diff --git a/src/RadhydroShockMultigroup/CMakeLists.txt b/src/RadhydroShockMultigroup/CMakeLists.txt new file mode 100644 index 000000000..8c379acce --- /dev/null +++ b/src/RadhydroShockMultigroup/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_radhydro_shock_multigroup test_radhydro_shock_multigroup.cpp ../fextract.cpp ../interpolate.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_shock_multigroup) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadhydroShockMultigroup COMMAND test_radhydro_shock_multigroup radshockMG.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp new file mode 100644 index 000000000..1640fb193 --- /dev/null +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -0,0 +1,410 @@ +/// \file test_radhydro_shock.cpp +/// \brief Defines a test problem for a multigroup gray-opacity non-equilibrium radiative shock. +/// + +#include + +#include "AMReX_Array.H" +#include "AMReX_BC_TYPES.H" + +#include "ArrayUtil.hpp" +#include "fextract.hpp" +#include "fundamental_constants.H" +#include "radiation_system.hpp" +#include "test_radhydro_shock_multigroup.hpp" +#include "valarray.hpp" +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif + +struct ShockProblem { +}; // dummy type to allow compile-type polymorphism via template specialization + +// parameters taken from Section 9.5 of Skinner et al. (2019) +// https://doi.org/10.3847/1538-4365/ab007f + +constexpr double a_rad = C::a_rad; // erg cm^-3 K^-4 +constexpr double c = C::c_light; // cm s^-1 +constexpr double k_B = C::k_B; // erg K^-1 + +constexpr double c_s0 = 1.73e7; // adiabatic sound speed [cm s^-1] +constexpr double kappa = 577.0; // "opacity" == rho*kappa [cm^-1] (!!) +constexpr double gamma_gas = (5. / 3.); +constexpr double c_v = k_B / ((C::m_p + C::m_e) * (gamma_gas - 1.0)); // specific heat [erg g^-1 K^-1] +constexpr double T0 = 2.18e6; // K +constexpr double rho0 = 5.69; // g cm^-3 +constexpr double v0 = 5.19e7; // cm s^-1 +constexpr double T1 = 7.98e6; // K [7.98297e6] +constexpr double rho1 = 17.1; // g cm^-3 [17.08233] +constexpr double v1 = 1.73e7; // cm s^-1 [1.72875e7] +constexpr double chat = 10.0 * (v0 + c_s0); // reduced speed of light + +constexpr double Erad0 = a_rad * (T0 * T0 * T0 * T0); // erg cm^-3 +constexpr double Egas0 = rho0 * c_v * T0; // erg cm^-3 +constexpr double Erad1 = a_rad * (T1 * T1 * T1 * T1); // erg cm^-3 +constexpr double Egas1 = rho1 * c_v * T1; // erg cm^-3 + +constexpr double shock_position = 0.0130; // 0.0132; // cm (shock position drifts to the right slightly during the simulation, so + // we initialize slightly to the left...) +constexpr double Lx = 0.01575; // cm + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = 8; +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = 0.; + static constexpr bool compute_v_over_c_terms = true; + static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz + static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{ + 1.00000000e+15, 3.16227766e+15, 1.00000000e+16, 3.16227766e+16, 1.00000000e+17, 3.16227766e+17, 1.00000000e+18, 3.16227766e+18, 1.00000000e+19}; +}; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = C::m_p + C::m_e; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = gamma_gas; +}; + +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) + -> quokka::valarray +{ + quokka::valarray kappaPVec{}; + for (int i = 0; i < nGroups_; ++i) { + kappaPVec[i] = kappa / rho; + } + return kappaPVec; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double /*Tgas*/) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, 0.0); +} + +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double +{ + return (1. / 3.); // Eddington approximation +} + +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 (!((bcr->lo(0) == amrex::BCType::ext_dir) || (bcr->hi(0) == amrex::BCType::ext_dir))) { // NOLINT + return; + } + +#if (AMREX_SPACEDIM == 1) + auto i = iv.toArray()[0]; + int const j = 0; + int const 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(); + amrex::GpuArray hi = box.hiVect3d(); + + auto const radBoundaries_g = RadSystem::radBoundaries_; + + if (i < lo[0]) { + const double Egas_L = Egas0 + 0.5 * rho0 * (v0 * v0); + const double xmom_L = consVar(lo[0], j, k, RadSystem::x1GasMomentum_index); + const double px_L = (xmom_L < (rho0 * v0)) ? xmom_L : (rho0 * v0); + + const auto Tgas0 = quokka::EOS::ComputeTgasFromEint(rho0, Egas_L); + auto radEnergyFractionsT0 = RadSystem::ComputePlanckEnergyFractions(radBoundaries_g, Tgas0); + + // x1 left side boundary -- constant + consVar(i, j, k, RadSystem::gasDensity_index) = rho0; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = 0.; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = px_L; // xmom_L; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas_L; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas_L - (px_L * px_L) / (2 * rho0); + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad0 * radEnergyFractionsT0[g]; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + } else if (i >= hi[0]) { + const double Egas_R = Egas1 + 0.5 * rho1 * (v1 * v1); + const double xmom_R = consVar(hi[0], j, k, RadSystem::x1GasMomentum_index); + const double px_R = (xmom_R > (rho1 * v1)) ? xmom_R : (rho1 * v1); + + const auto Tgas1 = quokka::EOS::ComputeTgasFromEint(rho1, Egas_R); + auto radEnergyFractionsT1 = RadSystem::ComputePlanckEnergyFractions(radBoundaries_g, Tgas1); + + // x1 right-side boundary -- constant + consVar(i, j, k, RadSystem::gasDensity_index) = rho1; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = 0.; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = px_R; // xmom_R; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas_R; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas_R - (px_R * px_R) / (2 * rho1); + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad1 * radEnergyFractionsT1[g]; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + } +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray const prob_lo = grid_elem.prob_lo_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + + // loop over the cell-centered quantities and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + + amrex::Real radEnergy = NAN; + amrex::Real x1RadFlux = NAN; + amrex::Real energy = NAN; + amrex::Real density = NAN; + amrex::Real x1Momentum = NAN; + amrex::Real temp = NAN; + quokka::valarray::nGroups> radEnergyFractions{}; + + if (x < shock_position) { + radEnergy = Erad0; + x1RadFlux = 0.0; + energy = Egas0 + 0.5 * rho0 * (v0 * v0); + density = rho0; + x1Momentum = rho0 * v0; + } else { + radEnergy = Erad1; + x1RadFlux = 0.0; + energy = Egas1 + 0.5 * rho1 * (v1 * v1); + density = rho1; + x1Momentum = rho1 * v1; + } + + temp = quokka::EOS::ComputeTgasFromEint(density, energy); + radEnergyFractions = RadSystem::ComputePlanckEnergyFractions(radBoundaries_g, temp); + + state_cc(i, j, k, RadSystem::gasDensity_index) = density; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = 0.; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = x1Momentum; + 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) = energy; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = energy - (x1Momentum * x1Momentum) / (2 * density); + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = radEnergy * radEnergyFractions[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = x1RadFlux; + 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; + } + }); +} + +auto problem_main() -> int +{ + // Problem parameters + const int max_timesteps = 2e4; + const double CFL_number = 0.4; + // const int nx = 512; + // const double initial_dtau = 1.0e-3; // dimensionless time + // const double max_dtau = 1.0e-3; // dimensionless time + // const double initial_dt = initial_dtau / c_s0; + // const double max_dt = max_dtau / c_s0; + const double max_time = 1.0e-9; // 9.08e-10; // s + + 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); // custom x1 + BCs_cc[n].setHi(0, amrex::BCType::ext_dir); // custom x1 + for (int i = 1; i < AMREX_SPACEDIM; ++i) { // x2- and x3- directions + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + sim.cflNumber_ = CFL_number; + sim.radiationCflNumber_ = CFL_number; + sim.maxTimesteps_ = max_timesteps; + sim.stopTime_ = max_time; + sim.plotfileInterval_ = -1; + + // run + sim.setInitialConditions(); + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); + int const nx = static_cast(position.size()); + + // Plot results + int status = 0; + if (amrex::ParallelDescriptor::IOProcessor()) { + std::vector xs(nx); + std::vector Trad(nx); + std::vector Tgas(nx); + std::vector Erad(nx); + std::vector Egas(nx); + + for (int i = 0; i < nx; ++i) { + const double x = Lx * ((i + 0.5) / static_cast(nx)); + xs.at(i) = x; // cm + + double Erad_t = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + Erad.at(i) = Erad_t / a_rad; // scaled + Trad.at(i) = std::pow(Erad_t / a_rad, 1. / 4.) / T0; // dimensionless + + const double Etot_t = values.at(RadSystem::gasEnergy_index)[i]; + const double rho = values.at(RadSystem::gasDensity_index)[i]; + const double x1GasMom = values.at(RadSystem::x1GasMomentum_index)[i]; + const double Ekin = (x1GasMom * x1GasMom) / (2.0 * rho); + + const double Egas_t = (Etot_t - Ekin); + Egas.at(i) = Egas_t; + Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho, Egas_t) / T0; // dimensionless + } + + // read in exact solution + + std::vector xs_exact; + std::vector Trad_exact; + std::vector Tmat_exact; + + std::string const filename = "../extern/LowrieEdwards/shock.txt"; + std::ifstream fstream(filename, std::ios::in); + + const double error_tol = 0.008; + double rel_error = NAN; + if (fstream.is_open()) { + + std::string header; + std::getline(fstream, header); + + for (std::string line; std::getline(fstream, line);) { + std::istringstream iss(line); + std::vector values; + + for (double value = NAN; iss >> value;) { + values.push_back(value); + } + auto x_val = values.at(0); // cm + auto Tmat_val = values.at(3); // dimensionless + auto Trad_val = values.at(4); // dimensionless + + if ((x_val > 0.0) && (x_val < Lx)) { + xs_exact.push_back(x_val); + Tmat_exact.push_back(Tmat_val); + Trad_exact.push_back(Trad_val); + } + } + + // compute error norm + std::vector Trad_interp(xs_exact.size()); + amrex::Print() << "xs min/max = " << xs[0] << ", " << xs[xs.size() - 1] << std::endl; + amrex::Print() << "xs_exact min/max = " << xs_exact[0] << ", " << xs_exact[xs_exact.size() - 1] << std::endl; + + interpolate_arrays(xs_exact.data(), Trad_interp.data(), static_cast(xs_exact.size()), xs.data(), Trad.data(), + static_cast(xs.size())); + + double err_norm = 0.; + double sol_norm = 0.; + for (size_t i = 0; i < xs_exact.size(); ++i) { + err_norm += std::abs(Trad_interp[i] - Trad_exact[i]); + sol_norm += std::abs(Trad_exact[i]); + } + + rel_error = err_norm / sol_norm; + amrex::Print() << "Error norm = " << err_norm << std::endl; + amrex::Print() << "Solution norm = " << sol_norm << std::endl; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + } + + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + +#ifdef HAVE_PYTHON + std::vector xs_scaled(xs.size()); + std::vector xs_exact_scaled(xs_exact.size()); + for (size_t i = 0; i < xs.size(); ++i) { + xs_scaled.at(i) = xs.at(i) / Lx; + } + for (size_t i = 0; i < xs_exact.size(); ++i) { + xs_exact_scaled.at(i) = xs_exact.at(i) / Lx; + } + + // plot results + int const s = 48; // stride + std::map Trad_args; + Trad_args["label"] = "radiation"; + Trad_args["color"] = "C1"; + matplotlibcpp::plot(xs_scaled, Trad, Trad_args); + + if (fstream.is_open()) { + std::unordered_map Trad_exact_args; + // Trad_exact_args["label"] = "Trad (diffusion ODE)"; + Trad_exact_args["color"] = "C1"; + Trad_exact_args["marker"] = "o"; + // Trad_exact_args["edgecolors"] = "k"; + matplotlibcpp::scatter(strided_vector_from(xs_exact_scaled, s), strided_vector_from(Trad_exact, s), 10.0, Trad_exact_args); + } + + std::map Tgas_args; + Tgas_args["label"] = "gas"; + Tgas_args["color"] = "C2"; + matplotlibcpp::plot(xs_scaled, Tgas, Tgas_args); + + if (fstream.is_open()) { + std::unordered_map Tgas_exact_args; + // Tgas_exact_args["label"] = "Tmat (diffusion ODE)"; + Tgas_exact_args["color"] = "C2"; + Tgas_exact_args["marker"] = "o"; + // Tgas_exact_args["edgecolors"] = "k"; + matplotlibcpp::scatter(strided_vector_from(xs_exact_scaled, s), strided_vector_from(Tmat_exact, s), 10.0, Tgas_exact_args); + } + + matplotlibcpp::xlabel("length x (dimensionless)"); + matplotlibcpp::ylabel("temperature (dimensionless)"); + matplotlibcpp::legend(); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./radshock_multigroup_temperature.pdf"); +#endif + } + + return status; +} diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.hpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.hpp new file mode 100644 index 000000000..2afd18eab --- /dev/null +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.hpp @@ -0,0 +1,28 @@ +#ifndef TEST_RADHYDRO_SHOCK_MULTIGROUP_HPP_ // NOLINT +#define TEST_RADHYDRO_SHOCK_MULTIGROUP_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_radhydro_shock_multigroup.hpp +/// \brief Defines a test problem for a radiative shock. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers +#include "RadhydroSimulation.hpp" +#include "hydro_system.hpp" +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// function definitions +auto problem_main() -> int; + +#endif // TEST_RADHYDRO_SHOCK_MULTIGROUP_HPP_ diff --git a/src/planck_integral.hpp b/src/planck_integral.hpp index 1b29da33f..05068dc9a 100644 --- a/src/planck_integral.hpp +++ b/src/planck_integral.hpp @@ -244,6 +244,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto integrate_planck_from_0_to_x(const if (logx < LOG_X_MIN) { // y = x * x * x / 3.0; // 1st order y = (-4 + x) * x + 8 * std::log((2 + x) / 2); // 2nd order + // Y_INTERP_MIN is the minimum value returned from interpolate_planck_integral. To ensure y is monotonic with respect to x: // AMREX_ASSERT(y <= Y_INTERP_MIN); if (y > Y_INTERP_MIN) { y = Y_INTERP_MIN; diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index ae2437e78..b688e7d72 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -825,6 +825,104 @@ AMREX_GPU_DEVICE void RadSystem::ComputeRadPressure(TARRAY &F, double S = std::max(0.1, std::sqrt(Tnormal)); } +template +template +AMREX_GPU_DEVICE void RadSystem::ComputeRadPressure(TARRAY &F, double &S, const double erad, const double Fx, const double Fy, const double Fz, + const double fx, const double fy, const double fz) +{ + // check that states are physically admissible + AMREX_ASSERT(erad > 0.0); + // AMREX_ASSERT(f < 1.0); // there is sometimes a small (<1%) flux + // limiting violation when using P1 AMREX_ASSERT(f_R < 1.0); + + auto f = std::sqrt(fx * fx + fy * fy + fz * fz); + std::array fvec = {fx, fy, fz}; + + // angle between interface and radiation flux \hat{n} + // If direction is undefined, just drop direction-dependent + // terms. + std::array n{}; + + for (int ii = 0; ii < 3; ++ii) { + n[ii] = (f > 0.) ? (fvec[ii] / f) : 0.; + } + + // compute radiation pressure tensors + const double chi = RadSystem::ComputeEddingtonFactor(f); + + AMREX_ASSERT((chi >= 1. / 3.) && (chi <= 1.0)); // NOLINT + + // diagonal term of Eddington tensor + const double Tdiag = (1.0 - chi) / 2.0; + + // anisotropic term of Eddington tensor (in the direction of the + // rad. flux) + const double Tf = (3.0 * chi - 1.0) / 2.0; + + // assemble Eddington tensor + std::array, 3> T{}; + std::array, 3> P{}; + + for (int ii = 0; ii < 3; ++ii) { + for (int jj = 0; jj < 3; ++jj) { + const double delta_ij = (ii == jj) ? 1 : 0; + T[ii][jj] = Tdiag * delta_ij + Tf * (n[ii] * n[jj]); + // compute the elements of the total radiation pressure + // tensor + P[ii][jj] = T[ii][jj] * erad; + } + } + + // frozen Eddington tensor approximation, following Balsara + // (1999) [JQSRT Vol. 61, No. 5, pp. 617–627, 1999], Eq. 46. + double Tnormal = NAN; + if constexpr (DIR == FluxDir::X1) { + Tnormal = T[0][0]; + } else if constexpr (DIR == FluxDir::X2) { + Tnormal = T[1][1]; + } else if constexpr (DIR == FluxDir::X3) { + Tnormal = T[2][2]; + } + + // compute fluxes F_L, F_R + // P_nx, P_ny, P_nz indicate components where 'n' is the direction of the + // face normal F_n is the radiation flux component in the direction of the + // face normal + double Fn = NAN; + double Pnx = NAN; + double Pny = NAN; + double Pnz = NAN; + + if constexpr (DIR == FluxDir::X1) { + Fn = Fx; + + Pnx = P[0][0]; + Pny = P[0][1]; + Pnz = P[0][2]; + } else if constexpr (DIR == FluxDir::X2) { + Fn = Fy; + + Pnx = P[1][0]; + Pny = P[1][1]; + Pnz = P[1][2]; + } else if constexpr (DIR == FluxDir::X3) { + Fn = Fz; + + Pnx = P[2][0]; + Pny = P[2][1]; + Pnz = P[2][2]; + } + + AMREX_ASSERT(Fn != NAN); + AMREX_ASSERT(Pnx != NAN); + AMREX_ASSERT(Pny != NAN); + AMREX_ASSERT(Pnz != NAN); + + F = {Fn, Pnx, Pny, Pnz}; + + S = std::max(0.1, std::sqrt(Tnormal)); +} + template template void RadSystem::ComputeFluxes(array_t &x1Flux_in, array_t &x1FluxDiffusive_in, amrex::Array4 const &x1LeftState_in, @@ -1317,6 +1415,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::GpuArray Frad_t0{}; amrex::GpuArray Frad_t1{}; amrex::GpuArray dMomentum{0., 0., 0.}; + std::array gasMtm{}; + double realFourPiB = NAN; T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); auto realFourPiB = c * ComputeThermalRadiation(T_gas, radBoundaries_g_copy); diff --git a/tests/RadhydroPulse.in b/tests/RadhydroPulse.in index 3f24bc52a..51ba0ecfb 100644 --- a/tests/RadhydroPulse.in +++ b/tests/RadhydroPulse.in @@ -16,6 +16,9 @@ amr.v = 0 # verbosity in Amr 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 +amr.max_grid_size_x = 64 +amr.max_grid_size_y = 4 +amr.max_grid_size_z = 4 do_reflux = 0 do_subcycle = 0 diff --git a/tests/radshockMG.in b/tests/radshockMG.in new file mode 100644 index 000000000..63ec0aeca --- /dev/null +++ b/tests/radshockMG.in @@ -0,0 +1,24 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 0.01575 0.001575 1.0 +geometry.is_periodic = 0 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 128 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this +amr.max_grid_size_x = 128 +amr.max_grid_size_y = 4 +amr.max_grid_size_z = 4 + +do_reflux = 0 +do_subcycle = 0