PopIII Simulation

src/CMakeLists.txt

@@ -106,6 +106,9 @@ set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SO
 #endif()
 #link_libraries(QuokkaObj)
 
+add_subdirectory(PopIII)
+
+
 add_subdirectory(StarCluster)
 add_subdirectory(PrimordialChem)
 
@@ -154,4 +157,3 @@ add_subdirectory(ShockCloud)
 add_subdirectory(PassiveScalar)
 add_subdirectory(FCQuantities)
 add_subdirectory(SphericalCollapse)
-add_subdirectory(NSCBC)

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!");

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_CURRENT_SOURCE_DIR}/perturbation.py --kmin=2 --kmax=${result} --size=${first_digits} --alpha=2 --f_solenoidal=1.0 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/TurbDataReader.cpp

@@ -0,0 +1,119 @@
+//==============================================================================
+// 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 TurbDataReader.cpp
+/// \brief Reads turbulent driving fields generated as cubic HDF5 arrays.
+///
+
+#include "TurbDataReader.hpp"
+#include "AMReX_Arena.H"
+#include "AMReX_BLassert.H"
+#include "AMReX_Print.H"
+#include "AMReX_TableData.H"
+#include <string>
+
+auto read_dataset(hid_t &file_id, char const *dataset_name) -> amrex::Table3D<double>
+{
+	// open dataset
+	hid_t dset_id = 0;
+	dset_id = H5Dopen2(file_id, dataset_name, H5P_DEFAULT);
+	AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dset_id != -1, "Can't open table!");
+
+	// get dimensions
+	hid_t const dspace = H5Dget_space(dset_id);
+	const int ndims = H5Sget_simple_extent_ndims(dspace);
+	std::vector<hsize_t> dims(ndims);
+	H5Sget_simple_extent_dims(dspace, dims.data(), nullptr);
+
+	uint64_t data_size = 1;
+	for (int idim = 0; idim < ndims; ++idim) {
+		data_size *= dims[idim];
+	}
+
+	// allocate array for dataset storage
+	auto *temp_data = new double[data_size];
+
+	// read dataset
+	herr_t status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data);
+	AMREX_ALWAYS_ASSERT_WITH_MESSAGE(status != -1, "Failed to read dataset!");
+
+	// close dataset
+	status = H5Dclose(dset_id);
+
+	// WARNING: Table3D uses column-major (Fortran-order) indexing, but HDF5
+	// tables use row-major (C-order) indexing!
+	amrex::GpuArray<int, 3> const lo{0, 0, 0};
+	amrex::GpuArray<int, 3> const hi{static_cast<int>(dims[0]), static_cast<int>(dims[1]), static_cast<int>(dims[2])};
+	auto table = amrex::Table3D<double>(temp_data, lo, hi);
+	return table;
+}
+
+void initialize_turbdata(turb_data &data, std::string &data_file)
+{
+	amrex::Print() << "Initializing turbulence data...\n";
+	amrex::Print() << fmt::format("data_file: {}.\n", data_file);
+
+	herr_t status = 0;
+	herr_t const h5_error = -1;
+
+	// open file
+	hid_t file_id = 0;
+	file_id = H5Fopen(data_file.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+	AMREX_ALWAYS_ASSERT_WITH_MESSAGE(file_id != h5_error, "Failed to open data file!");
+
+	data.dvx = read_dataset(file_id, "pertx");
+	data.dvy = read_dataset(file_id, "perty");
+	data.dvz = read_dataset(file_id, "pertz");
+
+	// close file
+	status = H5Fclose(file_id);
+}
+
+auto get_tabledata(amrex::Table3D<double> &in_t) -> amrex::TableData<double, 3>
+{
+	amrex::Array<int, 3> tlo{in_t.begin[0], in_t.begin[1], in_t.begin[2]};
+	amrex::Array<int, 3> thi{in_t.end[0] - 1, in_t.end[1] - 1, in_t.end[2] - 1};
+	amrex::TableData<double, 3> tableData(tlo, thi, amrex::The_Pinned_Arena());
+	auto h_table = tableData.table();
+
+	amrex::Print() << "Copying tableData on indices " << tlo << " to " << thi << ".\n";
+
+	// fill tableData
+	for (int i = tlo[0]; i <= thi[0]; ++i) {
+		for (int j = tlo[1]; j <= thi[1]; ++j) {
+			for (int k = tlo[2]; k <= thi[2]; ++k) {
+				h_table(i, j, k) = in_t(i, j, k);
+			}
+		}
+	}
+
+	return tableData;
+}
+
+auto computeRms(amrex::TableData<amrex::Real, 3> &dvx, amrex::TableData<amrex::Real, 3> &dvy, amrex::TableData<amrex::Real, 3> &dvz) -> amrex::Real
+{
+	amrex::Array<int, 3> tlo = dvx.lo();
+	amrex::Array<int, 3> thi = dvx.hi();
+	auto const &dvx_table = dvx.const_table();
+	auto const &dvy_table = dvy.const_table();
+	auto const &dvz_table = dvz.const_table();
+
+	// compute rms power
+	amrex::Real rms_sq = 0;
+	amrex::Long N = 0;
+	for (int i = tlo[0]; i <= thi[0]; ++i) {
+		for (int j = tlo[1]; j <= thi[1]; ++j) {
+			for (int k = tlo[2]; k <= thi[2]; ++k) {
+				amrex::Real const vx = dvx_table(i, j, k);
+				amrex::Real const vy = dvy_table(i, j, k);
+				amrex::Real const vz = dvz_table(i, j, k);
+				rms_sq += vx * vx + vy * vy + vz * vz;
+				++N;
+			}
+		}
+	}
+	rms_sq /= static_cast<amrex::Real>(N);
+	return std::sqrt(rms_sq);
+}

src/PopIII/TurbDataReader.hpp

@@ -0,0 +1,59 @@
+#ifndef TURBDATAREADER_HPP_ // NOLINT
+#define TURBDATAREADER_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 TurbDataReader.hpp
+/// \brief Reads turbulent driving fields generated as cubic HDF5 arrays.
+///
+
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <iostream>
+#include <string>
+
+#include "fmt/core.h"
+#include <H5Dpublic.h>
+#include <H5Ppublic.h>
+#include <hdf5.h>
+
+#include "AMReX.H"
+#include "AMReX_Array.H"
+#include "AMReX_BLassert.H"
+#include "AMReX_TableData.H"
+
+/* HDF5 definitions */
+
+#define HDF5_FILE_I4 H5T_STD_I32BE
+#define HDF5_FILE_I8 H5T_STD_I64BE
+#define HDF5_FILE_R4 H5T_IEEE_F32BE
+#define HDF5_FILE_R8 H5T_IEEE_F64BE
+#define HDF5_FILE_B8 H5T_STD_B8BE
+
+#define HDF5_I4 H5T_NATIVE_INT
+#define HDF5_I8 H5T_NATIVE_LLONG
+#define HDF5_R4 H5T_NATIVE_FLOAT
+#define HDF5_R8 H5T_NATIVE_DOUBLE
+#define HDF5_R16 H5T_NATIVE_LDOUBLE
+
+// Cooling table storage
+
+using turb_data = struct turb_data {
+	// values
+	amrex::Table3D<double> dvx;
+	amrex::Table3D<double> dvy;
+	amrex::Table3D<double> dvz;
+};
+
+void initialize_turbdata(turb_data &data, std::string &data_file);
+
+auto read_dataset(hid_t &file_id, char const *dataset_name) -> amrex::Table3D<double>;
+
+auto get_tabledata(amrex::Table3D<double> &in_t) -> amrex::TableData<double, 3>;
+
+auto computeRms(amrex::TableData<amrex::Real, 3> &dvx, amrex::TableData<amrex::Real, 3> &dvy, amrex::TableData<amrex::Real, 3> &dvz) -> amrex::Real;
+
+#endif // TURBDATAREADER_HPP_

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