From c7b83bc0da78e57572fa972eb9ffd449ef61e57c Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 25 Jan 2024 19:37:07 -0500 Subject: [PATCH] prototype --- src/AgoraGalaxy/CMakeLists.txt | 8 ++ src/AgoraGalaxy/galaxy.cpp | 166 +++++++++++++++++++++++++++++++++ src/AgoraGalaxy/galaxy.hpp | 22 +++++ src/CMakeLists.txt | 5 +- tests/AgoraGalaxy.in | 34 +++++++ 5 files changed, 233 insertions(+), 2 deletions(-) create mode 100644 src/AgoraGalaxy/CMakeLists.txt create mode 100644 src/AgoraGalaxy/galaxy.cpp create mode 100644 src/AgoraGalaxy/galaxy.hpp create mode 100644 tests/AgoraGalaxy.in diff --git a/src/AgoraGalaxy/CMakeLists.txt b/src/AgoraGalaxy/CMakeLists.txt new file mode 100644 index 000000000..cfb8e6bd9 --- /dev/null +++ b/src/AgoraGalaxy/CMakeLists.txt @@ -0,0 +1,8 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(agora_galaxy galaxy.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(agora_galaxy) + endif() + + add_test(NAME AgoraGalaxy COMMAND agora_galaxy AgoraGalaxy.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/AgoraGalaxy/galaxy.cpp b/src/AgoraGalaxy/galaxy.cpp new file mode 100644 index 000000000..74df7fc8f --- /dev/null +++ b/src/AgoraGalaxy/galaxy.cpp @@ -0,0 +1,166 @@ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2024 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file galaxy.cpp +/// \brief Defines a simulation using the AGORA isolated galaxy initial conditions. +/// + +#include + +#include "AMReX_BC_TYPES.H" +#include "AMReX_DistributionMapping.H" +#include "AMReX_Geometry.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" + +#include "EOS.hpp" +#include "RadhydroSimulation.hpp" +#include "fundamental_constants.H" +#include "galaxy.hpp" +#include "hydro_system.hpp" + +struct AgoraGalaxy { +}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = 5. / 3.; + static constexpr double mean_molecular_weight = C::m_u; + static constexpr double boltzmann_constant = C::k_B; +}; + +template <> struct HydroSystem_Traits { + static constexpr bool reconstruct_eint = true; +}; + +template <> struct Physics_Traits { + static constexpr bool is_hydro_enabled = true; + static constexpr bool is_radiation_enabled = false; + static constexpr bool is_mhd_enabled = false; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr int nGroups = 1; // number of radiation groups +}; + +template <> struct SimulationData { + std::vector radius{}; + std::vector vcirc{}; +}; + +template <> void RadhydroSimulation::preCalculateInitialConditions() +{ + // 1. read in circular velocity table + // 2. copy to GPU + // 3. save interpolator object + // + // TODO(bwibking): implement. +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::GpuArray dx = grid_elem.dx_; + const amrex::Array4 &state_cc = grid_elem.array_; + const amrex::GpuArray prob_lo = grid_elem.prob_lo_; + const amrex::Real gamma = quokka::EOS_Traits::gamma; + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + // Cartesian coordinates + 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]; + + // cylindrical coordinates + amrex::Real const R = std::sqrt(std::pow(x, 2) + std::pow(y, 2)); + amrex::Real const theta = std::atan2(x, y); + + // compute double exponential density profile + double const rho = 1.0e-22; // g cm^{-3} + + // interpolate circular velocity based on radius of cell center + double const vcirc = 0; + double const vx = vcirc * std::cos(theta); + double const vy = vcirc * std::sin(theta); + double const vz = 0; + double const vsq = vx * vx + vy * vy + vz * vz; + + // compute temperature + double T = NAN; + if (R < 20.0e3 * C::parsec) { + T = 1.0e4; // K + } else { + T = 1.0e6; // K + } + const double Eint = quokka::EOS::ComputeEintFromTgas(rho, T); + + state_cc(i, j, k, HydroSystem::density_index) = rho; + state_cc(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; + state_cc(i, j, k, HydroSystem::x2Momentum_index) = rho * vy; + state_cc(i, j, k, HydroSystem::x3Momentum_index) = rho * vz; + state_cc(i, j, k, HydroSystem::energy_index) = Eint + 0.5 * rho * vsq; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = Eint; + }); +} + +template <> void RadhydroSimulation::createInitialParticles() +{ + // read particles from ASCII file + const int nreal_extra = 4; // mass vx vy vz + CICParticles->SetVerbose(1); + CICParticles->InitFromAsciiFile("AgoraGalaxy_particles.txt", nreal_extra, nullptr); +} + +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 == "gpot") { + const int ncomp = ncomp_cc_in; + auto const &phi_arr = phi[lev].const_arrays(); + auto output = mf.arrays(); + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept { output[bx](i, j, k, ncomp) = phi_arr[bx](i, j, k); }); + } +} + +auto problem_main() -> int +{ + auto isNormalComp = [=](int n, int dim) { + if ((n == HydroSystem::x1Momentum_index) && (dim == 0)) { + return true; + } + if ((n == HydroSystem::x2Momentum_index) && (dim == 1)) { + return true; + } + if ((n == HydroSystem::x3Momentum_index) && (dim == 2)) { + return true; + } + return false; + }; + + 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) { + if (isNormalComp(n, i)) { + BCs_cc[n].setLo(i, amrex::BCType::reflect_odd); + BCs_cc[n].setHi(i, amrex::BCType::reflect_odd); + } else { + BCs_cc[n].setLo(i, amrex::BCType::reflect_even); + BCs_cc[n].setHi(i, amrex::BCType::reflect_even); + } + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + sim.doPoissonSolve_ = 1; // enable self-gravity + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + const int status = 0; + return status; +} diff --git a/src/AgoraGalaxy/galaxy.hpp b/src/AgoraGalaxy/galaxy.hpp new file mode 100644 index 000000000..7d67532f1 --- /dev/null +++ b/src/AgoraGalaxy/galaxy.hpp @@ -0,0 +1,22 @@ +#ifndef TEST_BINARY_ORBIT_HPP_ // NOLINT +#define TEST_BINARY_ORBIT_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 binary_orbit.hpp +/// \brief Defines a test problem for a binary orbit +/// + +// external headers +#include + +// internal headers +#include "hydro_system.hpp" +#include "interpolate.hpp" + +// function definitions +auto problem_main() -> int; + +#endif // TEST_BINARY_ORBIT_HPP_ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4d007bdd8..f066647de 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -if(QUOKKA_PYTHON) +if(QUOKKA_PYTHON) message(STATUS "Building Quokka with Python support (to disable, add '-DQUOKKA_PYTHON=OFF' to the CMake command-line options)") find_package(Python COMPONENTS Interpreter Development NumPy) else() @@ -115,7 +115,7 @@ link_libraries(yaml-cpp::yaml-cpp) include(CTest) #create an object library for files that should be compiled with all test problems -set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" +set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp" "${openPMDSources}" "${gamma_law_sources}") #we don't use it anymore because it gives issues on Cray systems @@ -176,3 +176,4 @@ add_subdirectory(NSCBC) add_subdirectory(StarCluster) add_subdirectory(PrimordialChem) add_subdirectory(BinaryOrbitCIC) +add_subdirectory(AgoraGalaxy) diff --git a/tests/AgoraGalaxy.in b/tests/AgoraGalaxy.in new file mode 100644 index 000000000..19bf49341 --- /dev/null +++ b/tests/AgoraGalaxy.in @@ -0,0 +1,34 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = -2.5e13 -2.5e13 -2.5e13 +geometry.prob_hi = 2.5e13 2.5e13 2.5e13 +geometry.is_periodic = 0 0 0 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 1 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 32 32 32 +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 = 32 # at least 128 for GPUs +amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells +amr.grid_eff = 1.0 + +cfl = 0.3 +max_timesteps = 8000 +stop_time = 5.0e7 # s + +do_reflux = 1 +do_subcycle = 0 +do_cic_particles = 1 # turns on CIC particles + +ascent_interval = -1 +checkpoint_interval = -1 +plottime_interval = 2.0e5 # s +derived_vars = gpot