Merge branch 'development' into chongchong/multigroup-adv-pulse

quokka-astro · Jan 7, 2024 · ff57787 · ff57787
2 parents 386c6da + 5451412
commit ff57787
Show file tree

Hide file tree

Showing 38 changed files with 790 additions and 17 deletions.
diff --git 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 → docs/tests/attach/hydro_shuosher.png b/docs/tests/hydro_shuosher.png → docs/tests/attach/hydro_shuosher.png
diff --git a/docs/tests/hydro_sms.png → docs/tests/attach/hydro_sms.png b/docs/tests/hydro_sms.png → docs/tests/attach/hydro_sms.png
diff --git a/docs/tests/radcoupling.png → docs/tests/attach/radcoupling.png b/docs/tests/radcoupling.png → docs/tests/attach/radcoupling.png
diff --git a/docs/tests/radcoupling_rsla.png → docs/tests/attach/radcoupling_rsla.png b/docs/tests/radcoupling_rsla.png → 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
diff --git a/docs/tests/attach/radhydro_pulse_density.png b/docs/tests/attach/radhydro_pulse_density.png
diff --git a/docs/tests/attach/radhydro_pulse_temperature-1.png b/docs/tests/attach/radhydro_pulse_temperature-1.png
diff --git a/docs/tests/attach/radhydro_pulse_temperature.png b/docs/tests/attach/radhydro_pulse_temperature.png
diff --git a/docs/tests/attach/radhydro_pulse_velocity-1.png b/docs/tests/attach/radhydro_pulse_velocity-1.png
diff --git a/docs/tests/attach/radhydro_pulse_velocity.png b/docs/tests/attach/radhydro_pulse_velocity.png
diff --git a/docs/tests/attach/radhydro_uniform_advecting_temperature-nobeta.png b/docs/tests/attach/radhydro_uniform_advecting_temperature-nobeta.png
diff --git a/docs/tests/attach/radhydro_uniform_advecting_temperature.png b/docs/tests/attach/radhydro_uniform_advecting_temperature.png
diff --git a/docs/tests/attach/radhydro_uniform_advecting_velocity-nobeta.png b/docs/tests/attach/radhydro_uniform_advecting_velocity-nobeta.png
diff --git a/docs/tests/attach/radhydro_uniform_advecting_velocity.png b/docs/tests/attach/radhydro_uniform_advecting_velocity.png
diff --git a/docs/tests/radshock_cgs_temperature.png → ...tests/attach/radshock_cgs_temperature.png b/docs/tests/radshock_cgs_temperature.png → ...tests/attach/radshock_cgs_temperature.png
diff --git 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.
diff --git 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
 
diff --git a/docs/tests/radhydro_pulse.rst 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
@@ -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
@@ -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.
diff --git 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.
diff --git 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.
diff --git a/extern/Microphysics b/extern/Microphysics
diff --git a/extern/amrex b/extern/amrex
diff --git a/extern/fmt b/extern/fmt
diff --git a/extern/openPMD-api b/extern/openPMD-api
diff --git a/extern/yaml-cpp b/extern/yaml-cpp
diff --git 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
@@ -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
@@ -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
@@ -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)