PopIII Simulation

src/CMakeLists.txt

@@ -125,6 +125,12 @@ set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SO
 #endif()
 #link_libraries(QuokkaObj)
 
+add_subdirectory(PopIII)
+
+
+add_subdirectory(StarCluster)
+add_subdirectory(PrimordialChem)
+
 add_subdirectory(Advection)
 add_subdirectory(Advection2D)
 add_subdirectory(AdvectionSemiellipse)
@@ -173,6 +179,4 @@ add_subdirectory(PassiveScalar)
 add_subdirectory(FCQuantities)
 add_subdirectory(SphericalCollapse)
 add_subdirectory(NSCBC)
-add_subdirectory(StarCluster)
-add_subdirectory(PrimordialChem)
 add_subdirectory(BinaryOrbitCIC)

src/Chemistry.hpp

@@ -31,9 +31,10 @@ namespace quokka::chemistry
 
 AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, const Real dt);
 
-template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed)
+template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed)
 {
 
+	BL_PROFILE("Chemistry::computeChemistry()");
 	for (amrex::MFIter iter(mf); iter.isValid(); ++iter) {
 		const amrex::Box &indexRange = iter.validbox();
 		auto const &state = mf.array(iter);
@@ -64,6 +65,11 @@ template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const R
 				chemstate.xn[nn] = inmfracs[nn];
 			}
 
+			// dont do chemistry in cells with densities below the minimum density specified
+			if (rho < min_density_allowed) {
+				return;
+			}
+
 			// stop the test if we have reached very high densities
 			if (rho > max_density_allowed) {
 				amrex::Abort("Density exceeded max_density_allowed!");
@@ -77,11 +83,19 @@ template <typename problem_t> void computeChemistry(amrex::MultiFab &mf, const R
 			// call the EOS to set initial internal energy e
 			eos(eos_input_re, chemstate);
 
+			if (dt < 0) {
+				amrex::Abort("Cannot do chemistry with dt < 0!");
+			}
+
 			// do the actual integration
 			// do it in .cpp so that it is not built at compile time for all tests
 			// which would otherwise slow down compilation due to the large RHS file
 			chemburner(chemstate, dt);
 
+			if (std::isnan(chemstate.xn[0]) || std::isnan(chemstate.rho)) {
+				amrex::Abort("Burner returned NAN");
+			}
+
 			if (!chemstate.success) {
 				amrex::Abort("VODE integration was unsuccessful!");
 			}

src/PopIII/CMakeLists.txt

@@ -0,0 +1,77 @@
+if (AMReX_SPACEDIM EQUAL 3)
+    # Define a custom target that runs the Python script to produce the input perturbations file
+
+
+    set(microphysics_network_name "primordial_chem") #this will override network_name to primordial_chem for this directory only
+    setup_target_for_microphysics_compilation(${microphysics_network_name} "${CMAKE_CURRENT_BINARY_DIR}/")
+
+    #use the BEFORE keyword so that these files get priority in compilation for targets in this directory
+    #this is critical to ensure the correct Microphysics files are linked to primordial chem target
+    include_directories(BEFORE ${primordial_chem_dirs} "${CMAKE_CURRENT_BINARY_DIR}/" "includes/extern_parameters.H" "includes/network_properties.H")
+
+    add_executable(popiii popiii.cpp ../TurbDataReader.cpp ../main.cpp ../GrackleDataReader.cpp ../CloudyCooling.cpp ../Chemistry.cpp ${primordial_chem_sources})
+    target_compile_definitions(popiii PUBLIC PRIMORDIAL_CHEM) #this will add #define PRIMORDIAL_CHEM
+
+    if(AMReX_GPU_BACKEND MATCHES "CUDA")
+        setup_target_for_cuda_compilation(popiii)
+    endif()
+
+    #need h5py to run perturbations.py file below
+    execute_process(
+        COMMAND Python3::Interpreter -c "h5py"
+        RESULT_VARIABLE EXIT_CODE
+        OUTPUT_QUIET
+    )
+
+    #re-run cmake if this file is changed
+    set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS ${CMAKE_SOURCE_DIR}/tests/PopIII.in)
+
+    # Read the 'PopIII.in' file line by line
+    file(READ ${CMAKE_SOURCE_DIR}/tests/PopIII.in pop_in_contents)
+
+    # Split the file contents into lines
+    string(REPLACE "\n" ";" pop_in_lines "${pop_in_contents}")
+
+    # Iterate through the lines and look for lines starting with 'amr.n_cell'
+    foreach(line IN LISTS pop_in_lines)
+        if (line MATCHES "^amr\\.n_cell.*")
+            set(match_line "${line}")
+            break()
+        endif()
+    endforeach()
+
+    if (match_line)
+        message(STATUS "picked line is ${match_line}")
+        # Use a regular expression to find consecutive digit
+        string(REGEX MATCHALL "[0-9]+" digits_list "${match_line}")
+        # Iterate through the matched digits and extract the first one
+        list(GET digits_list 0 first_digits)
+        # Check if matched digits were found
+        if (first_digits)
+           #first_digits give the box size, but you also want box size / 2 as kmax in the script below
+           math(EXPR int_first_digits "${first_digits}")
+           # Check if the conversion was successful
+           if (int_first_digits)
+               # Divide the integer by two
+               math(EXPR result "${int_first_digits} / 2")
+           else()
+               message(FATAL_ERROR "Conversion to integer failed, so cannot find kmax!")
+           endif()
+           message(STATUS "Will create initial perturbations of box size ${first_digits} and kmax ${result}")
+
+        else()
+           message(FATAL_ERROR "No box size found to create initial perturbations for!")
+        endif()
+    else()
+        message(FATAL_ERROR "No matching line found in ${CMAKE_SOURCE_DIR}/tests/PopIII.in!")
+    endif()
+
+
+    add_test(NAME ComputePerturbations COMMAND python3 ${CMAKE_SOURCE_DIR}/src/perturbation.py --kmin=2 --kmax=${result} --size=${first_digits} --alpha=1.8 --f_solenoidal=0.66667 WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    add_test(NAME PopIII COMMAND popiii PopIII.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    set_tests_properties(ComputePerturbations PROPERTIES FIXTURES_SETUP PopIII_fixture)
+    set_tests_properties(PopIII PROPERTIES FIXTURES_REQUIRED PopIII_fixture)
+
+    # AMR test only works on Setonix because Gadi and avatar do not have enough memory per GPU
+    # add_test(NAME PopIIIAMR COMMAND popiii popiii_AMR.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()

src/PopIII/ascent_actions.yaml

@@ -0,0 +1,30 @@
+- 
+  action: "add_pipelines"
+  pipelines: 
+    pl1:
+      f2:
+        type: "slice"
+        params: 
+          point: 
+            x: 0.
+            y: 0.
+            z: 0.
+          normal: 
+            x: 0.0
+            y: 0.0
+            z: 1.0
+- 
+  action: "add_scenes"
+  scenes: 
+    s1: 
+      plots: 
+        p1: 
+          type: "pseudocolor"
+          field: "log_density"
+          pipeline: "pl1"
+      renders:
+        r1:
+          image_prefix: "dens%05d"
+          annotations: "true"
+          camera:
+            zoom: 1.5