Merge branch 'development' into BenWibking/agora-galaxy

quokka-astro · Feb 3, 2024 · cc31803 · cc31803
2 parents 03e240d + 3455fc8
commit cc31803
Show file tree

Hide file tree

Showing 9 changed files with 662 additions and 74 deletions.
diff --git a/.ci/azure-pipelines-aarch64-debug.yml b/.ci/azure-pipelines-aarch64-debug.yml
@@ -30,7 +30,7 @@ jobs:
   - task: CMake@1
     displayName: 'Run CTest'
     inputs:
-      cmakeArgs: '-E chdir . ctest -j 2 -T Test --output-on-failure'
+      cmakeArgs: '-E chdir . ctest -j 2 -T Test -E "RadhydroPulseMG*" --output-on-failure'
 
   - task: PublishTestResults@2
     inputs:

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -172,6 +172,7 @@ add_subdirectory(RadhydroShockCGS)
 add_subdirectory(RadhydroShockMultigroup)
 add_subdirectory(RadhydroUniformAdvecting)
 add_subdirectory(RadhydroPulse)
+add_subdirectory(RadhydroPulseMG)
 
 add_subdirectory(ODEIntegration)
 add_subdirectory(Cooling)

diff --git a/src/RadhydroPulse/CMakeLists.txt b/src/RadhydroPulse/CMakeLists.txt
@@ -6,4 +6,4 @@ if (AMReX_SPACEDIM EQUAL 1)
   endif(AMReX_GPU_BACKEND MATCHES "CUDA")
 
   add_test(NAME RadhydroPulse COMMAND test_radhydro_pulse RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
-endif()
+endif()
diff --git a/src/RadhydroPulseMG/CMakeLists.txt b/src/RadhydroPulseMG/CMakeLists.txt
@@ -0,0 +1,9 @@
+if (AMReX_SPACEDIM EQUAL 1)
+  add_executable(test_radhydro_pulse_MG test_radhydro_pulse_MG.cpp ../fextract.cpp ${QuokkaObjSources})
+
+  if(AMReX_GPU_BACKEND MATCHES "CUDA")
+      setup_target_for_cuda_compilation(test_radhydro_pulse_MG)
+  endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+  add_test(NAME RadhydroPulseMG COMMAND test_radhydro_pulse_MG RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()
diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp
diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp b/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp
@@ -0,0 +1,21 @@
+#ifndef TEST_RADHYDRO_PULSE_MG_HPP_ // NOLINT
+#define TEST_RADHYDRO_PULSE_MG_HPP_
+/// \file test_radiation_marshak.hpp
+/// \brief Defines a test problem for radiation in the static diffusion regime.
+///
+
+// external headers
+#ifdef HAVE_PYTHON
+#include "matplotlibcpp.h"
+#endif
+#include <fmt/format.h>
+#include <fstream>
+
+// internal headers
+
+#include "interpolate.hpp"
+#include "radiation_system.hpp"
+
+// function definitions
+
+#endif // TEST_RADHYDRO_PULSE_MG_HPP_
diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp
@@ -40,6 +40,7 @@ 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 Erad_floor_ = Erad0 * 1e-12;
 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
@@ -63,7 +64,7 @@ template <> struct RadSystem_Traits<ShockProblem> {
 	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 double Erad_floor = Erad_floor_;
 	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<double, Physics_Traits<ShockProblem>::nGroups + 1> radBoundaries{

diff --git a/src/planck_integral.hpp b/src/planck_integral.hpp
@@ -245,8 +245,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto integrate_planck_from_0_to_x(const
 		// 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;
+		} else if (y < 0.) {
+			y = 0.;
 		}
 	} else if (logx >= LOG_X_MAX) {
 		return 1.0;