diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml index 3e2fd32e8..13994d366 100644 --- a/.github/workflows/clang-tidy.yml +++ b/.github/workflows/clang-tidy.yml @@ -23,7 +23,7 @@ jobs: config_file: src/.clang-tidy build_dir: build apt_packages: libopenmpi-dev,libhdf5-mpi-dev,python3-dev,python3-numpy,python3-matplotlib - cmake_command: cmake . -B build -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DQUOKKA_PYTHON=ON -DQUOKKA_OPENPMD=ON -DopenPMD_USE_ADIOS2=OFF + cmake_command: cmake . -B build -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DQUOKKA_PYTHON=ON -DQUOKKA_OPENPMD=ON -DopenPMD_USE_ADIOS2=OFF -DAMReX_SPACEDIM=3 split_workflow: true # Uploads an artefact containing clang_fixes.json diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7620f9726..036211169 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -24,4 +24,4 @@ ci: autoupdate_commit_msg: '[pre-commit.ci] pre-commit autoupdate' autoupdate_schedule: weekly skip: [] - submodules: false \ No newline at end of file + submodules: false diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 689728792..91e394b1f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -174,13 +174,14 @@ add_subdirectory(RadhydroUniformAdvecting) add_subdirectory(RadhydroPulse) add_subdirectory(RadhydroPulseMG) -add_subdirectory(ODEIntegration) +add_subdirectory(BinaryOrbitCIC) add_subdirectory(Cooling) -add_subdirectory(ShockCloud) -add_subdirectory(PassiveScalar) add_subdirectory(FCQuantities) -add_subdirectory(SphericalCollapse) add_subdirectory(NSCBC) -add_subdirectory(StarCluster) +add_subdirectory(ODEIntegration) +add_subdirectory(PassiveScalar) add_subdirectory(PrimordialChem) -add_subdirectory(BinaryOrbitCIC) +add_subdirectory(PopIII) +add_subdirectory(ShockCloud) +add_subdirectory(StarCluster) +add_subdirectory(SphericalCollapse) diff --git a/src/Chemistry.hpp b/src/Chemistry.hpp index 1a30a0c99..31a4248aa 100644 --- a/src/Chemistry.hpp +++ b/src/Chemistry.hpp @@ -31,13 +31,18 @@ namespace quokka::chemistry AMREX_GPU_DEVICE void chemburner(burn_t &chemstate, const Real dt); -template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed) +template void computeChemistry(amrex::MultiFab &mf, const Real dt, const Real max_density_allowed, const Real min_density_allowed) { + const BL_PROFILE("Chemistry::computeChemistry()"); for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { const amrex::Box &indexRange = iter.validbox(); auto const &state = mf.array(iter); + if (dt < 0) { + amrex::Abort("Cannot do chemistry with dt < 0!"); + } + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real rho = state(i, j, k, HydroSystem::density_index); const Real xmom = state(i, j, k, HydroSystem::x1Momentum_index); @@ -64,6 +69,11 @@ template 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!"); @@ -82,6 +92,10 @@ template void computeChemistry(amrex::MultiFab &mf, const R // 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!"); } diff --git a/src/PopIII/CMakeLists.txt b/src/PopIII/CMakeLists.txt new file mode 100644 index 000000000..b379db66f --- /dev/null +++ b/src/PopIII/CMakeLists.txt @@ -0,0 +1,78 @@ +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 "${openPMDSources}" ../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 + #purpose of this code is to check whether the h5py python package is installed + 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() diff --git a/src/PopIII/ascent_actions.yaml b/src/PopIII/ascent_actions.yaml new file mode 100644 index 000000000..6515c4b84 --- /dev/null +++ b/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 diff --git a/src/PopIII/popiii.cpp b/src/PopIII/popiii.cpp new file mode 100644 index 000000000..a4deb6ca2 --- /dev/null +++ b/src/PopIII/popiii.cpp @@ -0,0 +1,490 @@ +//============================================================================== +// 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 popiii.cpp +/// \brief Defines a test problem for Pop III star formation. +/// Author: Piyush Sharda (Leiden University, 2023) +/// +#include +#include +#include +#include +#include + +#include "AMReX.H" +#include "AMReX_Arena.H" +#include "AMReX_BC_TYPES.H" +#include "AMReX_BLProfiler.H" +#include "AMReX_BLassert.H" +#include "AMReX_Config.H" +#include "AMReX_FabArray.H" +#include "AMReX_FabArrayUtility.H" +#include "AMReX_GpuDevice.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParallelContext.H" +#include "AMReX_ParallelDescriptor.H" +#include "AMReX_ParmParse.H" +#include "AMReX_Print.H" +#include "AMReX_REAL.H" +#include "AMReX_SPACE.H" +#include "AMReX_TableData.H" + +#include "RadhydroSimulation.hpp" +#include "SimulationData.hpp" +#include "TurbDataReader.hpp" +#include "hydro_system.hpp" +#include "popiii.hpp" +#include "radiation_system.hpp" + +#include "actual_eos_data.H" +#include "burn_type.H" +#include "eos.H" +#include "extern_parameters.H" +#include "network.H" + +using amrex::Real; + +struct PopIII { +}; + +template <> struct HydroSystem_Traits { + static constexpr bool reconstruct_eint = true; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = NumSpec; // number of chemical species + static constexpr int numPassiveScalars = numMassScalars + 0; // we only have mass scalars + static constexpr bool is_radiation_enabled = false; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = 1; +}; + +template <> struct SimulationData { + // real-space perturbation fields + amrex::TableData dvx; + amrex::TableData dvy; + amrex::TableData dvz; + amrex::Real dv_rms_generated{}; + amrex::Real dv_rms_target{}; + amrex::Real rescale_factor{}; + + // cloud parameters + amrex::Real R_sphere{}; + amrex::Real numdens_init{}; + amrex::Real omega_sphere{}; + + amrex::Real small_temp{}; + amrex::Real small_dens{}; + amrex::Real temperature{}; + amrex::Real primary_species_1{}; + amrex::Real primary_species_2{}; + amrex::Real primary_species_3{}; + amrex::Real primary_species_4{}; + amrex::Real primary_species_5{}; + amrex::Real primary_species_6{}; + amrex::Real primary_species_7{}; + amrex::Real primary_species_8{}; + amrex::Real primary_species_9{}; + amrex::Real primary_species_10{}; + amrex::Real primary_species_11{}; + amrex::Real primary_species_12{}; + amrex::Real primary_species_13{}; + amrex::Real primary_species_14{}; +}; + +template <> void RadhydroSimulation::preCalculateInitialConditions() +{ + + // initialize microphysics routines + init_extern_parameters(); + + // parmparse species and temperature + amrex::ParmParse const pp("primordial_chem"); + userData_.small_temp = 1e1; + pp.query("small_temp", userData_.small_temp); + + userData_.small_dens = 1e-60; + pp.query("small_dens", userData_.small_dens); + + userData_.temperature = 1e1; + pp.query("temperature", userData_.temperature); + + userData_.primary_species_1 = 1.0e0_rt; + userData_.primary_species_2 = 0.0e0_rt; + userData_.primary_species_3 = 0.0e0_rt; + userData_.primary_species_4 = 0.0e0_rt; + userData_.primary_species_5 = 0.0e0_rt; + userData_.primary_species_6 = 0.0e0_rt; + userData_.primary_species_7 = 0.0e0_rt; + userData_.primary_species_8 = 0.0e0_rt; + userData_.primary_species_9 = 0.0e0_rt; + userData_.primary_species_10 = 0.0e0_rt; + userData_.primary_species_11 = 0.0e0_rt; + userData_.primary_species_12 = 0.0e0_rt; + userData_.primary_species_13 = 0.0e0_rt; + userData_.primary_species_14 = 0.0e0_rt; + + pp.query("primary_species_1", userData_.primary_species_1); + pp.query("primary_species_2", userData_.primary_species_2); + pp.query("primary_species_3", userData_.primary_species_3); + pp.query("primary_species_4", userData_.primary_species_4); + pp.query("primary_species_5", userData_.primary_species_5); + pp.query("primary_species_6", userData_.primary_species_6); + pp.query("primary_species_7", userData_.primary_species_7); + pp.query("primary_species_8", userData_.primary_species_8); + pp.query("primary_species_9", userData_.primary_species_9); + pp.query("primary_species_10", userData_.primary_species_10); + pp.query("primary_species_11", userData_.primary_species_11); + pp.query("primary_species_12", userData_.primary_species_12); + pp.query("primary_species_13", userData_.primary_species_13); + pp.query("primary_species_14", userData_.primary_species_14); + + eos_init(userData_.small_temp, userData_.small_dens); + network_init(); + + static bool isSamplingDone = false; + if (!isSamplingDone) { + // read perturbations from file + turb_data turbData; + amrex::ParmParse const pp("perturb"); + std::string turbdata_filename; + pp.query("filename", turbdata_filename); + initialize_turbdata(turbData, turbdata_filename); + + // copy to pinned memory + auto pinned_dvx = get_tabledata(turbData.dvx); + auto pinned_dvy = get_tabledata(turbData.dvy); + auto pinned_dvz = get_tabledata(turbData.dvz); + + // compute normalisation + userData_.dv_rms_generated = computeRms(pinned_dvx, pinned_dvy, pinned_dvz); + amrex::Print() << "rms dv = " << userData_.dv_rms_generated << "\n"; + + Real rms_dv_target = NAN; + pp.query("rms_velocity", rms_dv_target); + const Real rms_dv_actual = userData_.dv_rms_generated; + userData_.rescale_factor = rms_dv_target / rms_dv_actual; + + // copy to GPU + userData_.dvx.resize(pinned_dvx.lo(), pinned_dvx.hi()); + userData_.dvx.copy(pinned_dvx); + + userData_.dvy.resize(pinned_dvy.lo(), pinned_dvy.hi()); + userData_.dvy.copy(pinned_dvy); + + userData_.dvz.resize(pinned_dvz.lo(), pinned_dvz.hi()); + userData_.dvz.copy(pinned_dvz); + + isSamplingDone = true; + } +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // set initial conditions + amrex::GpuArray const dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + amrex::GpuArray prob_hi = grid_elem.prob_hi_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + std::array numdens = {-1.0}; + + for (int n = 1; n <= NumSpec; ++n) { + switch (n) { + + case 1: + numdens[n - 1] = userData_.primary_species_1; + break; + case 2: + numdens[n - 1] = userData_.primary_species_2; + break; + case 3: + numdens[n - 1] = userData_.primary_species_3; + break; + case 4: + numdens[n - 1] = userData_.primary_species_4; + break; + case 5: + numdens[n - 1] = userData_.primary_species_5; + break; + case 6: + numdens[n - 1] = userData_.primary_species_6; + break; + case 7: + numdens[n - 1] = userData_.primary_species_7; + break; + case 8: + numdens[n - 1] = userData_.primary_species_8; + break; + case 9: + numdens[n - 1] = userData_.primary_species_9; + break; + case 10: + numdens[n - 1] = userData_.primary_species_10; + break; + case 11: + numdens[n - 1] = userData_.primary_species_11; + break; + case 12: + numdens[n - 1] = userData_.primary_species_12; + break; + case 13: + numdens[n - 1] = userData_.primary_species_13; + break; + case 14: + numdens[n - 1] = userData_.primary_species_14; + break; + + default: + amrex::Abort("Cannot initialize number density for chem specie"); + break; + } + } + + amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); + amrex::Real const y0 = prob_lo[1] + 0.5 * (prob_hi[1] - prob_lo[1]); + amrex::Real const z0 = prob_lo[2] + 0.5 * (prob_hi[2] - prob_lo[2]); + + // cloud parameters + const double R_sphere = userData_.R_sphere; + const double omega_sphere = userData_.omega_sphere; + const double renorm_amp = userData_.rescale_factor; + const double numdens_init = userData_.numdens_init; + const double core_temp = userData_.temperature; + + auto const &dvx = userData_.dvx.const_table(); + auto const &dvy = userData_.dvy.const_table(); + auto const &dvz = userData_.dvz.const_table(); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; + amrex::Real const y = prob_lo[1] + (j + static_cast(0.5)) * dx[1]; + amrex::Real const z = prob_lo[2] + (k + static_cast(0.5)) * dx[2]; + amrex::Real const r = std::sqrt(std::pow(x - x0, 2) + std::pow(y - y0, 2) + std::pow(z - z0, 2)); + amrex::Real const distxy = std::sqrt(std::pow(x - x0, 2) + std::pow(y - y0, 2)); + + eos_t state; + amrex::Real rhotot = 0.0; + + for (int n = 0; n < NumSpec; ++n) { + state.xn[n] = numdens[n] * numdens_init; + rhotot += state.xn[n] * spmasses[n]; // spmasses contains the masses of all species, defined in EOS + } + + // normalize -- just in case + std::array mfracs = {-1.0}; + Real msum = 0.0; + for (int n = 0; n < NumSpec; ++n) { + mfracs[n] = state.xn[n] * spmasses[n] / rhotot; + msum += mfracs[n]; + } + + for (int n = 0; n < NumSpec; ++n) { + mfracs[n] /= msum; + // use the normalized mass fractions to obtain the corresponding number densities + state.xn[n] = mfracs[n] * rhotot / spmasses[n]; + } + + // amrex::Print() << "cell " << i << j << k << " " << rhotot << " " << numdens_init << " " << numdens[0] << std::endl; + + amrex::Real const phi = atan2((y - y0), (x - x0)); + + double vx = renorm_amp * dvx(i, j, k); + double vy = renorm_amp * dvy(i, j, k); + double const vz = renorm_amp * dvz(i, j, k); + + // calculate eos params for the core + state.rho = rhotot; + state.T = core_temp; + eos(eos_input_rt, state); + + if (r <= R_sphere) { + // add rotation to vx and vy + vx += (-1.0) * distxy * omega_sphere * std::sin(phi); + vy += distxy * omega_sphere * std::cos(phi); + + } else { + // re-calculate eos params outside the core, using pressure equilibrium (so, pressure within the core = pressure outside) + state.rho = 0.01 * rhotot; + state.p = state.p; + eos(eos_input_rp, state); + } + + // call the EOS to set initial internal energy e + amrex::Real const e = state.rho * state.e; + + // amrex::Print() << "cell " << i << j << k << " " << state.rho << " " << state.T << " " << e << std::endl; + + state_cc(i, j, k, HydroSystem::density_index) = state.rho; + state_cc(i, j, k, HydroSystem::x1Momentum_index) = state.rho * vx; + state_cc(i, j, k, HydroSystem::x2Momentum_index) = state.rho * vy; + state_cc(i, j, k, HydroSystem::x3Momentum_index) = state.rho * vz; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = e; + + Real const Egas = RadSystem::ComputeEgasFromEint(state.rho, state.rho * vx, state.rho * vy, state.rho * vz, e); + state_cc(i, j, k, HydroSystem::energy_index) = Egas; + + for (int nn = 0; nn < NumSpec; ++nn) { + state_cc(i, j, k, HydroSystem::scalar0_index + nn) = mfracs[nn] * state.rho; // we use partial densities and not mass fractions + } + }); +} + +template <> void RadhydroSimulation::ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real /*time*/, int /*ngrow*/) +{ + + // read-in jeans length refinement runtime params + amrex::ParmParse const pp("jeansRefine"); + int N_cells = 0; + pp.query("ncells", N_cells); // inverse of the 'Jeans number' [Truelove et al. (1997)] + Real jeans_density_threshold = NAN; + pp.query("density_threshold", jeans_density_threshold); + + const amrex::Real G = Gconst_; + const amrex::Real dx = geom[lev].CellSizeArray()[0]; + + auto const &prob_lo = geom[lev].ProbLoArray(); + auto const &prob_hi = geom[lev].ProbHiArray(); + + for (amrex::MFIter mfi(state_new_cc_[lev]); mfi.isValid(); ++mfi) { + const amrex::Box &box = mfi.validbox(); + const auto state = state_new_cc_[lev].const_array(mfi); + const auto tag = tags.array(mfi); + const int nidx = HydroSystem::density_index; + + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx; + amrex::Real const y = prob_lo[1] + (j + static_cast(0.5)) * dx; + amrex::Real const z = prob_lo[2] + (k + static_cast(0.5)) * dx; + + Real const rho = state(i, j, k, nidx); + Real const pressure = HydroSystem::ComputePressure(state, i, j, k); + amrex::GpuArray::numMassScalars> massScalars = RadSystem::ComputeMassScalars(state, i, j, k); + + amrex::Real const cs = quokka::EOS::ComputeSoundSpeed(rho, pressure, massScalars); + + const amrex::Real l_Jeans = cs * std::sqrt(M_PI / (G * rho)); + // add a density criterion for refinement so that no initial refinement is ever triggered outside the core + // typically, a density threshold ~ initial core density works well + if (l_Jeans < (N_cells * dx) && rho > jeans_density_threshold) { + tag(i, j, k) = amrex::TagBox::SET; + } + }); + } +} + +template <> void RadhydroSimulation::ComputeDerivedVar(int lev, std::string const &dname, amrex::MultiFab &mf, const int ncomp_cc_in) const +{ + // compute derived variables and save in 'mf' + if (dname == "temperature") { + const int ncomp = ncomp_cc_in; + auto const &state = state_new_cc_[lev].const_arrays(); + auto output = mf.arrays(); + + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { + Real const rho = state[bx](i, j, k, HydroSystem::density_index); + amrex::Real const Eint = state[bx](i, j, k, HydroSystem::internalEnergy_index); + + amrex::GpuArray::numMassScalars> massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); + + output[bx](i, j, k, ncomp) = quokka::EOS::ComputeTgasFromEint(rho, Eint, massScalars); + }); + } + + if (dname == "pressure") { + + const int ncomp = ncomp_cc_in; + auto const &state = state_new_cc_[lev].const_arrays(); + auto output = mf.arrays(); + + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { + amrex::Real const Pgas = HydroSystem::ComputePressure(state[bx], i, j, k); + output[bx](i, j, k, ncomp) = Pgas; + }); + } + + if (dname == "velx") { + + const int ncomp = ncomp_cc_in; + auto const &state = state_new_cc_[lev].const_arrays(); + auto output = mf.arrays(); + + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { + Real const rho = state[bx](i, j, k, HydroSystem::density_index); + Real const xmom = state[bx](i, j, k, HydroSystem::x1Momentum_index); + output[bx](i, j, k, ncomp) = xmom / rho; + }); + } + + if (dname == "sound_speed") { + + const int ncomp = ncomp_cc_in; + auto const &state = state_new_cc_[lev].const_arrays(); + auto output = mf.arrays(); + + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { + Real const rho = state[bx](i, j, k, HydroSystem::density_index); + Real const pressure = HydroSystem::ComputePressure(state[bx], i, j, k); + amrex::GpuArray::numMassScalars> massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); + + amrex::Real const cs = quokka::EOS::ComputeSoundSpeed(rho, pressure, massScalars); + output[bx](i, j, k, ncomp) = cs; + }); + } +} + +auto problem_main() -> int +{ + // read problem parameters + amrex::ParmParse const pp("perturb"); + + // cloud radius + Real R_sphere{}; + pp.query("cloud_radius", R_sphere); + + // cloud density + Real numdens_init{}; + pp.query("cloud_numdens", numdens_init); + + // cloud angular velocity + Real omega_sphere{}; + pp.query("cloud_omega", omega_sphere); + + // boundary conditions + const int ncomp_cc = Physics_Indices::nvarTotal_cc; + amrex::Vector BCs_cc(ncomp_cc); + for (int n = 0; n < ncomp_cc; ++n) { + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + BCs_cc[n].setLo(i, amrex::BCType::foextrap); + BCs_cc[n].setHi(i, amrex::BCType::foextrap); + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + sim.doPoissonSolve_ = 1; // enable self-gravity + + sim.tempFloor_ = 2.73 * (30.0 + 1.0); + // sim.speedCeiling_ = 3e6; + + sim.userData_.R_sphere = R_sphere; + sim.userData_.numdens_init = numdens_init; + sim.userData_.omega_sphere = omega_sphere; + + sim.initDt_ = 1e6; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + int const status = 0; + return status; +} diff --git a/src/PopIII/popiii.hpp b/src/PopIII/popiii.hpp new file mode 100644 index 000000000..9d3695bf3 --- /dev/null +++ b/src/PopIII/popiii.hpp @@ -0,0 +1,23 @@ +#ifndef POPIII_HPP_ // NOLINT +#define POPIII_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 star_cluster.hpp +/// \brief Defines a test problem for Pop III star formation. +/// Author: Piyush Sharda (Leiden University, 2023) +/// + +// external headers +#include + +// internal headers +#include "hydro_system.hpp" +#include "interpolate.hpp" + +// function definitions +auto problem_main() -> int; + +#endif // POPIII_HPP_ diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index ae1e5dee3..2793efd38 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -123,6 +123,8 @@ template class RadhydroSimulation : public AMRSimulation::max(); + Real min_density_allowed = std::numeric_limits::min(); + quokka::cooling::cloudy_tables cloudyTables_; std::string coolingTableFilename_{}; @@ -374,7 +376,8 @@ template void RadhydroSimulation::readParmParse( { amrex::ParmParse hpp("primordial_chem"); hpp.query("enabled", enableChemistry_); - hpp.query("max_density_allowed", max_density_allowed); + hpp.query("max_density_allowed", max_density_allowed); // chemistry is not accurate for densities > 3e-6 + hpp.query("min_density_allowed", min_density_allowed); // don't do chemistry in cells with densities below the minimum density specified } #endif @@ -514,7 +517,7 @@ void RadhydroSimulation::addStrangSplitSourcesWithBuiltin(amrex::Mult #ifdef PRIMORDIAL_CHEM if (enableChemistry_ == 1) { // compute chemistry - quokka::chemistry::computeChemistry(state, dt, max_density_allowed); + quokka::chemistry::computeChemistry(state, dt, max_density_allowed, min_density_allowed); } #endif @@ -1053,7 +1056,7 @@ auto RadhydroSimulation::advanceHydroAtLevel(amrex::MultiFab &state_o #else // write AMReX plotfile // WriteSingleLevelPlotfile(CustomPlotFileName("debug_stage1_filled_state_old", istep[lev]+1), - // state_old_cc_tmp, componentNames_cc_, geom[lev], time, istep[lev]+1); + // state_old_cc_tmp, componentNames_cc_, geom[lev], time, istep[lev]+1); #endif } @@ -1568,7 +1571,7 @@ void RadhydroSimulation::subcycleRadiationAtLevel(int lev, amrex::Rea for (int i = 0; i < nsubSteps; ++i) { if (i > 0) { // since we are starting a new substep, we need to copy radiation state from - // new state vector to old state vector + // new state vector to old state vector // (this is not necessary for the i=0 substep because we have already swapped // the full hydro+radiation state vectors at the beginning of the level advance) swapRadiationState(state_old_cc_[lev], state_new_cc_[lev]); diff --git a/src/StarCluster/CMakeLists.txt b/src/StarCluster/CMakeLists.txt index 277525261..796f1f1c8 100644 --- a/src/StarCluster/CMakeLists.txt +++ b/src/StarCluster/CMakeLists.txt @@ -1,7 +1,7 @@ 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}) + add_executable(star_cluster star_cluster.cpp ../TurbDataReader.cpp ${QuokkaObjSources}) if(AMReX_GPU_BACKEND MATCHES "CUDA") setup_target_for_cuda_compilation(star_cluster) endif() @@ -12,7 +12,7 @@ if (AMReX_SPACEDIM EQUAL 3) 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 ComputeStarClusterPerturbations COMMAND python3 ${CMAKE_SOURCE_DIR}/src/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) diff --git a/src/StarCluster/TurbDataReader.cpp b/src/TurbDataReader.cpp similarity index 96% rename from src/StarCluster/TurbDataReader.cpp rename to src/TurbDataReader.cpp index 6834c17a1..cb5300b68 100644 --- a/src/StarCluster/TurbDataReader.cpp +++ b/src/TurbDataReader.cpp @@ -40,7 +40,7 @@ auto read_dataset(hid_t &file_id, char const *dataset_name) -> amrex::Table3D &in_t) -> amrex::TableData @@ -78,7 +78,7 @@ auto get_tabledata(amrex::Table3D &in_t) -> amrex::TableData amrex::TableData tableData(tlo, thi, amrex::The_Pinned_Arena()); auto h_table = tableData.table(); - amrex::Print() << "Copying tableData on indices " << tlo << " to " << thi << ".\n"; + // amrex::Print() << "Copying tableData on indices " << tlo << " to " << thi << ".\n"; // fill tableData for (int i = tlo[0]; i <= thi[0]; ++i) { diff --git a/src/StarCluster/TurbDataReader.hpp b/src/TurbDataReader.hpp similarity index 100% rename from src/StarCluster/TurbDataReader.hpp rename to src/TurbDataReader.hpp diff --git a/src/StarCluster/perturbation.py b/src/perturbation.py similarity index 99% rename from src/StarCluster/perturbation.py rename to src/perturbation.py index 888d1756c..b9a0f2717 100644 --- a/src/StarCluster/perturbation.py +++ b/src/perturbation.py @@ -262,7 +262,7 @@ def plot_spectrum1D(pertx, perty, pertz): print("You must choose f_solenoidal. See --help.") sys.exit(0) alpha = options.alpha -if alpha==None: +if alpha is None: print("You must choose a power law slope, alpha. See --help.") sys.exit(0) alpha = float(options.alpha) @@ -273,7 +273,7 @@ def plot_spectrum1D(pertx, perty, pertz): # data precision dtype = np.float64 # ratio of solenoidal to compressive components -if options.f_solenoidal=="None" or options.f_solenoidal==None: +if options.f_solenoidal is None or options.f_solenoidal is None: f_solenoidal = None else: f_solenoidal = min(max(float(options.f_solenoidal), 0.), 1.) diff --git a/tests/PopIII.in b/tests/PopIII.in new file mode 100644 index 000000000..3810ada3b --- /dev/null +++ b/tests/PopIII.in @@ -0,0 +1,108 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = -3.703e18 -3.703e18 -3.703e18 +geometry.prob_hi = 3.703e18 3.703e18 3.703e18 +geometry.is_periodic = 0 0 0 +#restartfile = "chk15360" +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 1 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 64 64 64 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 32 # grid size must be divisible by this +amr.max_grid_size = 128 # at least 128 for GPUs +amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells +amr.grid_eff = 0.7 # default + +hydro.reconstruction_order = 3 # PPM +cfl = 0.15 +max_timesteps = 10 +stop_time = 1e16 + +do_reflux = 1 +do_subcycle = 0 + +#ascent_interval = 50 +plotfile_interval = 50 #100 +checkpoint_interval = 200 + +perturb.cloud_radius = 3.086e18 #initial cloud radius in cm +perturb.cloud_omega = 2.016008E-14 #initial cloud angular velocity in s^-1 +perturb.cloud_numdens = 0.90861183E+004 #initial cloud number density in cm^-3 +perturb.rms_velocity = 1.8050e5 #initial cloud rms velocity (to drive turbulence) + +# params for jeans refinement +jeansRefine.ncells = 64 #refine on these many cells per jeans length +jeansRefine.jeans_density_threshold = 2e-20 #do not refine if density is less than this density + +# density floor for popiii +density_floor = 1e-25 + +# in quokka/src/StarCluster, generate with 'python3 perturbation.py --kmin=2 --kmax=64 --size=128 --alpha=2 --f_solenoidal=1.0' +# and put it in quokka/tests/ +perturb.filename = "zdrv.hdf5" + +derived_vars = temperature velx pressure sound_speed + +amrex.throw_exception = 0 +amrex.signal_handling = 1 + +primordial_chem.enabled = 1 +primordial_chem.temperature = 0.26415744E+003 + +primordial_chem.small_temp = 1.e1 +primordial_chem.small_dens = 1.e-60 +primordial_chem.max_density_allowed = 3e-6 +primordial_chem.min_density_allowed = 1e-21 +#format in krome fort22_wD.dat file: E H- D- H HE H2 HD D H+ HE+ H2+ D+ HD+ HE++ +#format in quokka: E H+ H H- D+ D H2+ D- H2 HD+ HD HE++ HE+ HE +primordial_chem.primary_species_1 = 0.88499253E-006 +primordial_chem.primary_species_2 = 0.88498062E-006 +primordial_chem.primary_species_3 = 0.99932238E+000 +primordial_chem.primary_species_4 = 0.54719550E-013 +primordial_chem.primary_species_5 = 0.21957612E-010 +primordial_chem.primary_species_6 = 0.29920413E-004 +primordial_chem.primary_species_7 = 0.58304958E-015 +primordial_chem.primary_species_8 = 0.22122496E-017 +primordial_chem.primary_species_9 = 0.38932607E-003 +primordial_chem.primary_species_10 = 0.36774691E-019 +primordial_chem.primary_species_11 = 0.79574711E-007 +primordial_chem.primary_species_12 = 0.39651766E-050 +primordial_chem.primary_species_13 = 0.24136647E-043 +primordial_chem.primary_species_14 = 0.77500001E-001 + + +# integrator runtime parameters +# are we using primordial chemistry? then we use number densities +integrator.use_number_densities = 1 +# we do not want to subtract the internal energy +integrator.subtract_internal_energy = 0 +# we do not want to clip species between 0 and 1 +integrator.do_species_clip = 0 +# minimum positive value of number densities +integrator.SMALL_X_SAFE = 1e-100 +integrator.burner_verbose = 0 + +# do you want to use the jacobian calculated in a previous step? +integrator.use_jacobian_caching = 0 +# integration will fail if the number density > X_reject_buffer*atol +integrator.X_reject_buffer = 1e100 +# Set which jacobian to use +# 1 = analytic jacobian +# 2 = numerical jacobian +integrator.jacobian = 1 + +# do you want to normalize abundances within VODE? (you don't!) +integrator.renormalize_abundances = 0 +# tolerances +integrator.rtol_spec = 1.0e-4 +integrator.atol_spec = 1.0e-4 + +#assumed redshift for Pop III star formation +network.redshift = 30.0