Starcluster example problem (#342)

* adjust multigrid parameters for gravity solve * add prototype star cluster example * reduce reltol for poisson solve * initialize with gaussian random field * add Jeans refinement * update params * change debugging settings * move settings to runtime params * add derived var * gpu fixes * gpu fix * fix OpenBCSolver init * add StarCluster AMR params file * delete old file * refactor code * read perturbations from file * fix TableData indices * increase max_retries * fix build on Crusher * add TinyProfiler regions * rename file content * define problem_main * generate zdrv.hdf5 * update file * update file * update microp * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fix spellings * reorder tests * reorder tests * run python script for starcluster test * add filename * add custom command to run python script * revert changes * use python to generate input file * copy input file to binary * reset maxtries to 4 * make Gconst a runtime parameter with a default value * copy zdrv.hdf5 to tests/ * remove StaarCluster AMR test from repo * remove sim.Gconst * updated microp * produce input file at runtime * check for h5py * synced microp; * update microp * clang-tidy comments * clang-tidy * fix bug * increase simulation time time to ~tff * set Gconst to unity --------- Co-authored-by: Ben Wibking <[email protected]> Co-authored-by: Ben Wibking <[email protected]> Co-authored-by: Piyush Sharda <[email protected]> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Piyush Sharda <[email protected]> Co-authored-by: Piyush Sharda <[email protected]>
quokka-astro · Aug 9, 2023 · fb96683 · fb96683
1 parent add2204
commit fb96683
Show file tree

Hide file tree

Showing 17 changed files with 934 additions and 14 deletions.
diff --git a/.gitmodules b/.gitmodules
@@ -10,4 +10,4 @@
 [submodule "extern/Microphysics"]
 	path = extern/Microphysics
 	url = https://github.com/psharda/Microphysics
-	branch = compute-cs
+	branch = development
diff --git a/cmake/crusher.cmake b/cmake/crusher.cmake
@@ -12,8 +12,7 @@ set(CMAKE_CXX_COMPILER "CC" CACHE PATH "")
 set(AMReX_GPU_BACKEND HIP CACHE STRING "")
 set(AMReX_AMD_ARCH gfx90a CACHE STRING "") # MI250X
 set(AMREX_GPUS_PER_NODE 8 CACHE STRING "")
-set(ENV{Ascent_DIR} "/ccs/home/wibking/ascent/ascent-v0.9.0-pre/install")
-set(ENV{Kokkos_DIR} "/ccs/home/wibking/ascent/kokkos-3.6.01/install/lib64/cmake/Kokkos")
+set(ENV{Ascent_DIR} "/ccs/home/wibking/ascent_crusher/install/ascent-develop")
 set(AMReX_CONDUIT ON BOOL STRING "")
 set(AMReX_ASCENT ON BOOL STRING "")
 option(QUOKKA_PYTHON OFF)
diff --git a/extern/Microphysics b/extern/Microphysics
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -105,7 +105,7 @@ set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SO
 #endif()
 #link_libraries(QuokkaObj)
 
-add_subdirectory(HydroShocktubeCMA)
+add_subdirectory(StarCluster)
 add_subdirectory(PrimordialChem)
 
 add_subdirectory(Advection)
@@ -120,6 +120,7 @@ add_subdirectory(HydroKelvinHelmholz)
 add_subdirectory(HydroLeblanc)
 add_subdirectory(HydroRichtmeyerMeshkov)
 add_subdirectory(HydroShocktube)
+add_subdirectory(HydroShocktubeCMA)
 add_subdirectory(HydroShuOsher)
 add_subdirectory(HydroSMS)
 add_subdirectory(HydroVacuum)
@@ -152,4 +153,3 @@ add_subdirectory(ShockCloud)
 add_subdirectory(PassiveScalar)
 add_subdirectory(FCQuantities)
 add_subdirectory(SphericalCollapse)
-
diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp
@@ -123,6 +123,8 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 
 	amrex::Long radiationCellUpdates_ = 0; // total number of radiation cell-updates
 
+	amrex::Real Gconst_ = C::Gconst; // gravitational constant G
+
 	// member functions
 	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc, amrex::Vector<amrex::BCRec> &BCs_fc) : AMRSimulation<problem_t>(BCs_cc, BCs_fc)
 	{
@@ -314,6 +316,12 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::readParmParse(
 		hpp.query("artificial_viscosity_coefficient", artificialViscosityK_);
 	}
 
+	// set gravity runtime parameter
+	{
+		amrex::ParmParse hpp("gravity");
+		hpp.query("Gconst", Gconst_);
+	}
+
 	// set cooling runtime parameters
 	{
 		amrex::ParmParse hpp("cooling");
@@ -606,10 +614,11 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::fillPoissonRhs
 	// NOTE: in the future, this should also deposit particle mass
 	auto const &state = state_new_cc_[lev].const_arrays();
 	auto rhs = rhs_mf.arrays();
+	const Real G = Gconst_;
 
 	amrex::ParallelFor(rhs_mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
 		// copy density to rhs_mf
-		rhs[bx](i, j, k) = 4.0 * M_PI * state[bx](i, j, k, HydroSystem<problem_t>::density_index);
+		rhs[bx](i, j, k) = 4.0 * M_PI * G * state[bx](i, j, k, HydroSystem<problem_t>::density_index);
 	});
 }
 
@@ -778,6 +787,7 @@ template <typename problem_t>
 void RadhydroSimulation<problem_t>::advanceHydroAtLevelWithRetries(int lev, amrex::Real time, amrex::Real dt_lev, amrex::YAFluxRegister *fr_as_crse,
 								   amrex::YAFluxRegister *fr_as_fine)
 {
+	BL_PROFILE_REGION("HydroSolver");
 	// timestep retries
 	const int max_retries = 4;
 	bool success = false;

diff --git a/src/SphericalCollapse/spherical_collapse.hpp b/src/SphericalCollapse/spherical_collapse.hpp
@@ -1,12 +1,12 @@
-#ifndef TEST_HYDRO3D_BLAST_HPP_ // NOLINT
-#define TEST_HYDRO3D_BLAST_HPP_
+#ifndef TEST_SPHERICAL_COLLAPSE_HPP_ // NOLINT
+#define TEST_SPHERICAL_COLLAPSE_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_hydro3d_blast.hpp
-/// \brief Defines a test problem for a 3D explosion.
+/// \file spherical_collapse.hpp
+/// \brief Defines a test problem for a pressureless spherical collapse
 ///
 
 // external headers
@@ -17,6 +17,6 @@
 #include "interpolate.hpp"
 
 // function definitions
-auto testproblem_hydro_sedov() -> int;
+auto problem_main() -> int;
 
-#endif // TEST_HYDRO3D_BLAST_HPP_
+#endif // TEST_SPHERICAL_COLLAPSE_HPP_
diff --git a/src/StarCluster/CMakeLists.txt b/src/StarCluster/CMakeLists.txt
@@ -0,0 +1,22 @@
+if (AMReX_SPACEDIM EQUAL 3)
+    # Define a custom target that runs the Python script to produce the input perturbations file
+
+    add_executable(star_cluster star_cluster.cpp TurbDataReader.cpp ${QuokkaObjSources})
+    if(AMReX_GPU_BACKEND MATCHES "CUDA")
+        setup_target_for_cuda_compilation(star_cluster)
+    endif()
+
+    execute_process(
+        COMMAND Python3::Interpreter -c "h5py"
+        RESULT_VARIABLE EXIT_CODE
+        OUTPUT_QUIET
+    )
+
+    add_test(NAME ComputeStarClusterPerturbations COMMAND python3 ${CMAKE_CURRENT_SOURCE_DIR}/perturbation.py --kmin=2 --kmax=64 --size=128 --alpha=2 --f_solenoidal=1.0 WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    add_test(NAME StarCluster COMMAND star_cluster StarCluster.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    set_tests_properties(ComputeStarClusterPerturbations PROPERTIES FIXTURES_SETUP    StarCluster_fixture)
+    set_tests_properties(StarCluster PROPERTIES FIXTURES_REQUIRED StarCluster_fixture)
+
+    # AMR test only works on Setonix because Gadi and avatar do not have enough memory per GPU
+    # add_test(NAME StarClusterAMR COMMAND star_cluster StarCluster_AMR.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()
diff --git a/src/StarCluster/TurbDataReader.cpp b/src/StarCluster/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);
+}
diff --git a/src/StarCluster/TurbDataReader.hpp b/src/StarCluster/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_
diff --git a/src/StarCluster/ascent_actions.yaml b/src/StarCluster/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