diff --git a/src/AdvectionSimulation.hpp b/src/AdvectionSimulation.hpp index dcb8959b2..9c5014bbd 100644 --- a/src/AdvectionSimulation.hpp +++ b/src/AdvectionSimulation.hpp @@ -72,6 +72,7 @@ template class AdvectionSimulation : public AMRSimulation

&initSumCons) override; @@ -155,6 +156,13 @@ template void AdvectionSimulation::setInitialCon // note: an implementation is only required if face-centered vars are used } +template void AdvectionSimulation::createInitialParticles() +{ + // default empty implementation + // user should implement using problem-specific template specialization + // note: an implementation is only required if particles are used +} + template void AdvectionSimulation::computeAfterTimestep() { // do nothing -- user should implement using problem-specific template specialization diff --git a/src/BinaryOrbitCIC/CMakeLists.txt b/src/BinaryOrbitCIC/CMakeLists.txt new file mode 100644 index 000000000..166171f73 --- /dev/null +++ b/src/BinaryOrbitCIC/CMakeLists.txt @@ -0,0 +1,8 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(binary_orbit binary_orbit.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(binary_orbit) + endif() + + add_test(NAME BinaryOrbitCIC COMMAND binary_orbit BinaryOrbit.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/BinaryOrbitCIC/binary_orbit.cpp b/src/BinaryOrbitCIC/binary_orbit.cpp new file mode 100644 index 000000000..b5782945e --- /dev/null +++ b/src/BinaryOrbitCIC/binary_orbit.cpp @@ -0,0 +1,197 @@ +//============================================================================== +// 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.cpp +/// \brief Defines a test problem for a binary orbit. +/// + +#include "AMReX.H" +#include "AMReX_Array.H" +#include "AMReX_BC_TYPES.H" +#include "AMReX_BLassert.H" +#include "AMReX_Config.H" +#include "AMReX_DistributionMapping.H" +#include "AMReX_FabArrayUtility.H" +#include "AMReX_Geometry.H" +#include "AMReX_GpuContainers.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParallelDescriptor.H" +#include "AMReX_ParmParse.H" +#include "AMReX_Print.H" + +#include "AMReX_REAL.H" +#include "AMReX_ccse-mpi.H" +#include "RadhydroSimulation.hpp" +#include "binary_orbit.hpp" +#include "hydro_system.hpp" +#include + +struct BinaryOrbit { +}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = 1.0; // isothermal + static constexpr double cs_isothermal = 1.3e7; // cm s^{-1} + 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 = false; +}; + +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 time{}; + std::vector dist{}; +}; + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + double const rho = 1.0e-22; // g cm^{-3} + state_cc(i, j, k, HydroSystem::density_index) = rho; + state_cc(i, j, k, HydroSystem::x1Momentum_index) = 0; + state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0; + state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0; + state_cc(i, j, k, HydroSystem::energy_index) = 0; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = 0; + }); +} + +template <> void RadhydroSimulation::createInitialParticles() +{ + // read particles from ASCII file + const int nreal_extra = 4; // mass vx vy vz + CICParticles->SetVerbose(1); + CICParticles->InitFromAsciiFile("BinaryOrbit_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); }); + } +} + +template <> void RadhydroSimulation::computeAfterTimestep() +{ + // every N cycles, save particle statistics + static int cycle = 1; + if (cycle % 10 == 0) { + // create single-box particle container + amrex::ParticleContainer analysisPC{}; + amrex::Box const box(amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, amrex::IntVect{AMREX_D_DECL(1, 1, 1)}); + amrex::Geometry const geom(box); + amrex::BoxArray const boxArray(box); + amrex::DistributionMapping const dmap(boxArray, 1); + analysisPC.Define(geom, dmap, boxArray); + analysisPC.copyParticles(*CICParticles); + // do we need to redistribute?? + + if (amrex::ParallelDescriptor::IOProcessor()) { + quokka::CICParticleIterator const pIter(analysisPC, 0); + if (pIter.isValid()) { // this returns false when there is more than 1 MPI rank (?) + amrex::Print() << "Computing particle statistics...\n"; + const amrex::Long np = pIter.numParticles(); + auto &particles = pIter.GetArrayOfStructs(); + + // copy particles from device to host + quokka::CICParticleContainer::ParticleType *pData = particles().data(); + amrex::Vector pData_h(np); + amrex::Gpu::copy(amrex::Gpu::deviceToHost, pData, pData + np, pData_h.begin()); // NOLINT + + // compute orbital elements + quokka::CICParticleContainer::ParticleType &p1 = pData_h[0]; + quokka::CICParticleContainer::ParticleType &p2 = pData_h[1]; + const amrex::ParticleReal dx = p1.pos(0) - p2.pos(0); + const amrex::ParticleReal dy = p1.pos(1) - p2.pos(1); + const amrex::ParticleReal dz = p1.pos(2) - p2.pos(2); + const amrex::ParticleReal dist = std::sqrt(dx * dx + dy * dy + dz * dz); + const amrex::ParticleReal dist0 = 6.25e12; // cm + const amrex::Real cell_dx0 = this->geom[0].CellSize(0); + + // save statistics + userData_.time.push_back(tNew_[0]); + userData_.dist.push_back((dist - dist0) / cell_dx0); + } + } + } + ++cycle; +} + +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 + sim.initDt_ = 1.0e3; // s + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // check max abs particle distance + double max_err = NAN; + if (amrex::ParallelDescriptor::IOProcessor() && (!sim.userData_.dist.empty())) { + auto result = std::max_element(sim.userData_.dist.begin(), sim.userData_.dist.end(), + [](amrex::ParticleReal a, amrex::ParticleReal b) { return std::abs(a) < std::abs(b); }); + max_err = std::abs(*result); + } + amrex::ParallelDescriptor::Bcast(&max_err, 1, amrex::ParallelDescriptor::Mpi_typemap::type(), amrex::ParallelDescriptor::ioProcessor, + amrex::ParallelDescriptor::Communicator()); + amrex::Print() << "max particle separation = " << max_err << " cell widths.\n"; + + int status = 1; + const double max_err_tol = 0.18; // max error tol in cell widths + if (max_err < max_err_tol) { + status = 0; + } + return status; +} diff --git a/src/BinaryOrbitCIC/binary_orbit.hpp b/src/BinaryOrbitCIC/binary_orbit.hpp new file mode 100644 index 000000000..7d67532f1 --- /dev/null +++ b/src/BinaryOrbitCIC/binary_orbit.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/BinaryOrbitCIC/particle_orbit_plot.py b/src/BinaryOrbitCIC/particle_orbit_plot.py new file mode 100644 index 000000000..d486b3bd7 --- /dev/null +++ b/src/BinaryOrbitCIC/particle_orbit_plot.py @@ -0,0 +1,42 @@ +# note: the AMReX frontend is broken in yt 4.3.0 +import yt +from yt.frontends.boxlib.data_structures import AMReXDataset + +def particle_dist(plotfiles): + t_arr = [] + err_arr = [] + d0 = 2.0 * 3.125e12 + + for pltfile in plotfiles: + ds = AMReXDataset(pltfile) + Lx = ds.domain_right_edge[0] - ds.domain_left_edge[0] + Nx = ds.domain_dimensions[0] + cell_dx = Lx/Nx + ad = ds.all_data() + x = ad["CIC_particles", "particle_position_x"] + y = ad["CIC_particles", "particle_position_y"] + z = ad["CIC_particles", "particle_position_z"] + dx = x[0] - x[1] + dy = y[0] - y[1] + dz = z[0] - z[1] + from math import sqrt + d = sqrt(dx*dx + dy*dy + dz*dz) + #fractional_err = (d-d0)/d0 + grid_err = (d-d0)/cell_dx + t_arr.append(float(ds.current_time) / 3.15e7) + err_arr.append(grid_err) + + return t_arr, err_arr + +import glob +files = glob.glob("plt*") +files = sorted(files) +t, err = particle_dist(files) + +import matplotlib.pyplot as plt +plt.figure(figsize=(6,4)) +plt.plot(t, err) +plt.xlabel("time (yr)") +plt.ylabel(r"$(d-d_0)/\Delta x$") +plt.tight_layout() +plt.savefig("orbit.png", dpi=150) diff --git a/src/CICParticles.hpp b/src/CICParticles.hpp new file mode 100644 index 000000000..6e4eaa500 --- /dev/null +++ b/src/CICParticles.hpp @@ -0,0 +1,48 @@ +#ifndef CICPARTICLES_HPP_ // NOLINT +#define CICPARTICLES_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 CICParticles.hpp +/// \brief Implements the particle container for gravitationally-interacting particles +/// + +#include "AMReX.H" +#include "AMReX_AmrParticles.H" +#include "AMReX_MultiFab.H" +#include "AMReX_MultiFabUtil.H" +#include "AMReX_ParIter.H" +#include "AMReX_ParticleInterpolators.H" +#include "AMReX_Particles.H" +namespace quokka +{ + +enum ParticleDataIdx { ParticleMassIdx = 0, ParticleVxIdx, ParticleVyIdx, ParticleVzIdx }; +constexpr int CICParticleRealComps = 4; // mass vx vy vz + +using CICParticleContainer = amrex::AmrParticleContainer; +using CICParticleIterator = amrex::ParIter; + +struct CICDeposition { + amrex::Real Gconst{}; + int start_part_comp{}; + int start_mesh_comp{}; + int num_comp{}; + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void operator()(const CICParticleContainer::ParticleType &p, amrex::Array4 const &rho, + amrex::GpuArray const &plo, + amrex::GpuArray const &dxi) const noexcept + { + amrex::ParticleInterpolator::Linear interp(p, plo, dxi); + interp.ParticleToMesh(p, rho, start_part_comp, start_mesh_comp, num_comp, + [=] AMREX_GPU_DEVICE(const CICParticleContainer::ParticleType &part, int comp) { + return 4.0 * M_PI * Gconst * part.rdata(comp); // weight by 4 pi G + }); + } +}; + +} // namespace quokka + +#endif // CICPARTICLES_HPP_ \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 590dd089d..4d007bdd8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -125,14 +125,10 @@ set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SO #endif() #link_libraries(QuokkaObj) -add_subdirectory(StarCluster) -add_subdirectory(PrimordialChem) - add_subdirectory(Advection) add_subdirectory(Advection2D) add_subdirectory(AdvectionSemiellipse) - add_subdirectory(HydroBlast2D) add_subdirectory(HydroBlast3D) add_subdirectory(HydroContact) @@ -177,3 +173,6 @@ add_subdirectory(PassiveScalar) add_subdirectory(FCQuantities) add_subdirectory(SphericalCollapse) add_subdirectory(NSCBC) +add_subdirectory(StarCluster) +add_subdirectory(PrimordialChem) +add_subdirectory(BinaryOrbitCIC) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 150c37cb4..ae1e5dee3 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -35,6 +35,7 @@ #include "AMReX_FabFactory.H" #include "AMReX_Geometry.H" #include "AMReX_GpuControl.H" +#include "AMReX_GpuDevice.H" #include "AMReX_GpuQualifiers.H" #include "AMReX_IArrayBox.H" #include "AMReX_IndexType.H" @@ -110,6 +111,7 @@ template class RadhydroSimulation : public AMRSimulation::WriteCheckpointFile; using AMRSimulation::GetData; using AMRSimulation::FillPatchWithData; + using AMRSimulation::Gconst_; using AMRSimulation::densityFloor_; using AMRSimulation::tempFloor_; @@ -146,8 +148,6 @@ template class RadhydroSimulation : public AMRSimulation &BCs_cc, amrex::Vector &BCs_fc) : AMRSimulation(BCs_cc, BCs_fc) { @@ -182,6 +182,7 @@ template class RadhydroSimulation : public AMRSimulation void RadhydroSimulation::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"); @@ -485,6 +480,13 @@ template void RadhydroSimulation::setInitialCond // note: an implementation is only required if face-centered vars are used } +template void RadhydroSimulation::createInitialParticles() +{ + // default empty implementation + // user should implement using problem-specific template specialization + // note: an implementation is only required if particles are used +} + template void RadhydroSimulation::computeAfterTimestep() { // do nothing -- user should implement if desired @@ -681,15 +683,16 @@ template void RadhydroSimulation::advanceSingleT template void RadhydroSimulation::fillPoissonRhsAtLevel(amrex::MultiFab &rhs_mf, const int lev) { // add hydro density to Poisson rhs - // 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 * G * state[bx](i, j, k, HydroSystem::density_index); + // *add* density to rhs_mf + // (N.B. particles **will not work** if you overwrite the density here!) + rhs[bx](i, j, k) += 4.0 * M_PI * G * state[bx](i, j, k, HydroSystem::density_index); }); + amrex::Gpu::streamSynchronizeAll(); } template void RadhydroSimulation::applyPoissonGravityAtLevel(amrex::MultiFab const &phi_mf, const int lev, const amrex::Real dt) diff --git a/src/SphericalCollapse/spherical_collapse.cpp b/src/SphericalCollapse/spherical_collapse.cpp index 1aeaf6967..c96d610d6 100644 --- a/src/SphericalCollapse/spherical_collapse.cpp +++ b/src/SphericalCollapse/spherical_collapse.cpp @@ -85,6 +85,19 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid }); } +template <> void RadhydroSimulation::createInitialParticles() +{ + // add particles at random positions in the box + const bool generate_on_root_rank = true; + const int iseed = 42; + const int num_particles = 1000; + const double total_particle_mass = 0.5; // about 0.1 of the total fluid mass + const double particle_mass = total_particle_mass / static_cast(num_particles); + + const quokka::CICParticleContainer::ParticleInitData pdata = {{particle_mass, 0, 0, 0}, {}, {}, {}}; // {mass vx vy vz}, empty, empty, empty + CICParticles->InitRandom(num_particles, iseed, pdata, generate_on_root_rank); +} + template <> void RadhydroSimulation::ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real /*time*/, int /*ngrow*/) { // tag cells for refinement @@ -110,12 +123,9 @@ template <> void RadhydroSimulation::ComputeDerivedVar(int lev, // compute derived variables and save in 'mf' if (dname == "gpot") { const int ncomp = ncomp_cc_in; - auto const &state = state_new_cc_[lev].const_arrays(); - auto const &phi_data_ref = phi[lev].const_arrays(); + 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_data_ref[bx](i, j, k, ncomp); }); + 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); }); } } diff --git a/src/simulation.hpp b/src/simulation.hpp index 085ac8516..14ad4d1d7 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -12,7 +12,9 @@ // c++ headers #include "AMReX_Interpolater.H" #include +#include #include +#include #include #if __has_include() #include @@ -63,6 +65,7 @@ namespace filesystem = experimental::filesystem; #include "AMReX_Vector.H" #include "AMReX_VisMF.H" #include "AMReX_YAFluxRegister.H" +#include "fundamental_constants.H" #include "physics_numVars.hpp" #include #include @@ -74,6 +77,7 @@ namespace filesystem = experimental::filesystem; #include #ifdef AMREX_PARTICLES +#include "CICParticles.hpp" #include #include #endif @@ -104,6 +108,8 @@ using namespace conduit; using namespace ascent; #endif +enum class ParticleStep { BeforePoissonSolve, AfterPoissonSolve }; + using variant_t = std::variant; namespace YAML @@ -229,6 +235,7 @@ template class AMRSimulation : public amrex::AmrCore virtual void preCalculateInitialConditions() = 0; virtual void setInitialConditionsOnGrid(quokka::grid grid_elem) = 0; virtual void setInitialConditionsOnGridFaceVars(quokka::grid grid_elem) = 0; + virtual void createInitialParticles() = 0; virtual void computeAfterTimestep() = 0; virtual void computeAfterEvolve(amrex::Vector &initSumCons) = 0; virtual void fillPoissonRhsAtLevel(amrex::MultiFab &rhs, int lev) = 0; @@ -329,6 +336,10 @@ template class AMRSimulation : public amrex::AmrCore [[nodiscard]] auto getOldMF_fc() const -> amrex::Vector> const &; [[nodiscard]] auto getNewMF_fc() const -> amrex::Vector> const &; + // particle functions + void kickParticlesAllLevels(amrex::Real dt); + void driftParticlesAllLevels(amrex::Real dt); + #ifdef AMREX_USE_ASCENT void AscentCustomActions(conduit::Node const &blueprintMesh); void RenderAscent(); @@ -378,11 +389,17 @@ template class AMRSimulation : public amrex::AmrCore amrex::Long cellUpdates_ = 0; amrex::Vector cellUpdatesEachLevel_; + // gravity + amrex::Real Gconst_ = C::Gconst; // gravitational constant G + // tracer particles #ifdef AMREX_PARTICLES - void InitParticles(); // create tracer particles + void InitParticles(); // create tracer particles + void InitCICParticles(); // create CIC particles int do_tracers = 0; + int do_cic_particles = 0; std::unique_ptr TracerPC; + std::unique_ptr CICParticles; #endif // external objects @@ -571,6 +588,9 @@ template void AMRSimulation::readParameters() // Default do_tracers = 0 (turns on/off tracer particles) pp.query("do_tracers", do_tracers); + // Default do_cic_particles = 0 (turns on/off CIC particles) + pp.query("do_cic_particles", do_cic_particles); + // Default suppress_output = 0 pp.query("suppress_output", suppress_output); @@ -608,6 +628,12 @@ template void AMRSimulation::readParameters() maxWalltime_ = 3600 * hours + 60 * minutes + seconds; amrex::Print() << fmt::format("Setting walltime limit to {} hours, {} minutes, {} seconds.\n", hours, minutes, seconds); } + + // set gravity runtime parameters + { + const amrex::ParmParse hpp("gravity"); + hpp.query("Gconst", Gconst_); + } } template void AMRSimulation::setInitialConditions() @@ -624,6 +650,9 @@ template void AMRSimulation::setInitialCondition if (do_tracers != 0) { InitParticles(); } + if (do_cic_particles != 0) { + InitCICParticles(); + } #endif if (checkpointInterval_ > 0) { @@ -830,14 +859,25 @@ template void AMRSimulation::evolve() amrex::ParallelDescriptor::Barrier(); // synchronize all MPI ranks computeTimestep(); + // do particle leapfrog (first kick at time t) + kickParticlesAllLevels(dt_[0]); + // hyperbolic advance over all levels + // (N.B. when AMR is enabled, regridding may happen during this function!) int lev = 0; // coarsest level const int iteration = 1; // this is the first call to advance level 'lev' timeStepWithSubcycling(lev, cur_time, iteration); + // drift particles from t to (t + dt) + // N.B.: MUST be done *before* Poisson solve at new time! + driftParticlesAllLevels(dt_[0]); + // elliptic solve over entire AMR grid (post-timestep) ellipticSolveAllLevels(dt_[0]); + // do particle leapfrog (second kick at t + dt) + kickParticlesAllLevels(dt_[0]); + cur_time += dt_[0]; ++cycleCount_; computeAfterTimestep(); @@ -981,7 +1021,20 @@ template void AMRSimulation::calculateGpotAllLev rhs[lev].define(grids[lev], dmap[lev], ncomp, nghost); phi[lev].setVal(0); // set initial guess to zero rhs[lev].setVal(0); + } + +#ifdef AMREX_PARTICLES + if (do_cic_particles != 0) { + // deposit particles using amrex::ParticleToMesh + amrex::ParticleToMesh(*CICParticles, amrex::GetVecOfPtrs(rhs), 0, finest_level, + quokka::CICDeposition{Gconst_, quokka::ParticleMassIdx, 0, 1}); + } +#endif + + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_ALWAYS_ASSERT(!rhs[lev].contains_nan()); fillPoissonRhsAtLevel(rhs[lev], lev); + AMREX_ALWAYS_ASSERT(!rhs[lev].contains_nan()); rhs_min = std::min(rhs_min, rhs[lev].min(0)); } @@ -990,6 +1043,11 @@ template void AMRSimulation::calculateGpotAllLev if (verbose) { amrex::Print() << "\n"; } + + // check for NaN + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_ALWAYS_ASSERT(!phi[lev].contains_nan()); // this fails when max_level=2 for SphericalCollapse + } } #endif } @@ -1021,6 +1079,126 @@ template void AMRSimulation::ellipticSolveAllLev #endif } +struct setFunctorParticleAccel { + AMREX_GPU_DEVICE void operator()(const amrex::IntVect &iv, amrex::Array4 const &dest, const int &dcomp, const int &numcomp, + amrex::GeometryData const &geom, const amrex::Real &time, const amrex::BCRec *bcr, int bcomp, + const int &orig_comp) const + { + amrex::ignore_unused(iv, dest, dcomp, numcomp, geom, time, bcr, bcomp, orig_comp); + } +}; + +template void AMRSimulation::kickParticlesAllLevels(const amrex::Real dt) +{ + // kick particles (do: vel[i] += 0.5 * dt * accel[i]) + + if (do_cic_particles != 0) { + // gravitational acceleration multifabs + amrex::Vector accel(finest_level + 1); + + // self-gravity in Quokka requires open boundary conditions, + // so we extrapolate the gravitational accelerations at physical boundaries + amrex::Vector accelBC(AMREX_SPACEDIM); + for (int j = 0; j < AMREX_SPACEDIM; ++j) { + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + accelBC[j].setLo(i, amrex::BCType::foextrap); + accelBC[j].setHi(i, amrex::BCType::foextrap); + } + } + + for (int lev = 0; lev <= finest_level; ++lev) { + // compute accelerations + accel[lev].define(boxArray(lev), DistributionMap(lev), AMREX_SPACEDIM, 1); + accel[lev].setVal(0.); + auto accel_arr = accel[lev].arrays(); + const auto &phi_arr = phi[lev].const_arrays(); + const auto dx_inv = geom[lev].InvCellSizeArray(); + const amrex::IntVect ng(0); + + // check for NaN + AMREX_ALWAYS_ASSERT(!phi[lev].contains_nan()); + + amrex::ParallelFor(accel[lev], ng, AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k, int n) { + // compute cell-centered acceleration -grad(phi) + if (n == 0) { + accel_arr[bx](i, j, k, n) = -0.5 * dx_inv[0] * (phi_arr[bx](i + 1, j, k) - phi_arr[bx](i - 1, j, k)); + } + if (n == 1) { + accel_arr[bx](i, j, k, n) = -0.5 * dx_inv[1] * (phi_arr[bx](i, j + 1, k) - phi_arr[bx](i, j - 1, k)); + } + if (n == 2) { + accel_arr[bx](i, j, k, n) = -0.5 * dx_inv[2] * (phi_arr[bx](i, j, k + 1) - phi_arr[bx](i, j, k - 1)); + } + }); + amrex::Gpu::streamSynchronizeAll(); + + // fill ghost cells for accel[lev] + amrex::GpuBndryFuncFab boundaryFunctor(setFunctorParticleAccel{}); + amrex::PhysBCFunct> fineBdryFunct(geom[lev], accelBC, boundaryFunctor); + + if (lev == 0) { + accel[lev].FillBoundary(geom[lev].periodicity()); + fineBdryFunct(accel[lev], 0, accel[lev].nComp(), accel[lev].nGrowVect(), 0., 0); + } else { + amrex::PhysBCFunct> coarseBdryFunct(geom[lev - 1], accelBC, boundaryFunctor); + amrex::InterpFromCoarseLevel(accel[lev], 0., accel[lev - 1], 0, 0, AMREX_SPACEDIM, geom[lev - 1], geom[lev], coarseBdryFunct, 0, + fineBdryFunct, 0, refRatio(lev - 1), getAmrInterpolaterCellCentered(), accelBC, 0); + } + + // check for NaN + AMREX_ALWAYS_ASSERT(!accel[lev].contains_nan(0, AMREX_SPACEDIM)); + AMREX_ALWAYS_ASSERT(!accel[lev].contains_nan()); + + // loop over boxes of particles on this level + for (quokka::CICParticleIterator pIter(*CICParticles, lev); pIter.isValid(); ++pIter) { + auto &particles = pIter.GetArrayOfStructs(); + quokka::CICParticleContainer::ParticleType *pData = particles().data(); + const amrex::Long np = pIter.numParticles(); + + amrex::Array4 const &accel_arr = accel[lev].array(pIter); + const auto plo = geom[lev].ProbLoArray(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int64_t idx) { + quokka::CICParticleContainer::ParticleType &p = pData[idx]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic) + amrex::ParticleInterpolator::Linear interp(p, plo, dx_inv); + interp.MeshToParticle( + p, accel_arr, 0, quokka::ParticleVxIdx, AMREX_SPACEDIM, + [=] AMREX_GPU_DEVICE(amrex::Array4 const &acc, int i, int j, int k, int comp) { + return acc(i, j, k, comp); // no weighting + }, + [=] AMREX_GPU_DEVICE(quokka::CICParticleContainer::ParticleType & p, int comp, amrex::Real acc_comp) { + // kick particle by updating its velocity + p.rdata(comp) += 0.5 * dt * static_cast(acc_comp); + }); + }); + } + } + } +} + +template void AMRSimulation::driftParticlesAllLevels(const amrex::Real dt) +{ + // drift all particles (do: pos[i] += dt * vel[i]) + + if (do_cic_particles != 0) { + for (int lev = 0; lev <= finest_level; ++lev) { + for (quokka::CICParticleIterator pIter(*CICParticles, lev); pIter.isValid(); ++pIter) { + auto &particles = pIter.GetArrayOfStructs(); + quokka::CICParticleContainer::ParticleType *pData = particles().data(); + const amrex::Long np = pIter.numParticles(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int64_t idx) { + quokka::CICParticleContainer::ParticleType &p = pData[idx]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic) + // update particle position + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + p.pos(i) += dt * p.rdata(quokka::ParticleVxIdx + i); + } + }); + } + } + } +} + // N.B.: This function actually works for subcycled or not subcycled, as long as // nsubsteps[lev] is set correctly. template void AMRSimulation::timeStepWithSubcycling(int lev, amrex::Real time, int iteration) @@ -1062,6 +1240,9 @@ template void AMRSimulation::timeStepWithSubcycl if (do_tracers != 0) { TracerPC->Redistribute(lev); } + if (do_cic_particles != 0) { + CICParticles->Redistribute(lev); + } #endif // do fix-up on all levels that have been re-gridded @@ -1118,7 +1299,7 @@ template void AMRSimulation::timeStepWithSubcycl } #ifdef AMREX_PARTICLES - // redistribute particles + // redistribute tracer particles if (do_tracers != 0) { int redistribute_ngrow = 0; if ((iteration < nsubsteps[lev]) || (lev == 0)) { @@ -1130,6 +1311,18 @@ template void AMRSimulation::timeStepWithSubcycl TracerPC->Redistribute(lev, TracerPC->finestLevel(), redistribute_ngrow); } } + // redistribute CIC particles + if (do_cic_particles != 0) { + int redistribute_ngrow = 0; + if ((iteration < nsubsteps[lev]) || (lev == 0)) { + if (lev == 0) { + redistribute_ngrow = 0; + } else { + redistribute_ngrow = iteration; + } + CICParticles->Redistribute(lev, CICParticles->finestLevel(), redistribute_ngrow); + } + } #endif } @@ -1747,6 +1940,18 @@ template void AMRSimulation::InitParticles() TracerPC->Redistribute(); } } + +template void AMRSimulation::InitCICParticles() +{ + if (do_cic_particles != 0) { + AMREX_ASSERT(CICParticles == nullptr); + CICParticles = std::make_unique(this); + + CICParticles->SetVerbose(0); + createInitialParticles(); + CICParticles->Redistribute(); + } +} #endif // get plotfile name @@ -1970,7 +2175,10 @@ template void AMRSimulation::WritePlotFile() #ifdef AMREX_PARTICLES // write particles if (do_tracers != 0) { - TracerPC->WritePlotFile(plotfilename, "particles"); + TracerPC->WritePlotFile(plotfilename, "tracer_particles"); + } + if (do_cic_particles != 0) { + CICParticles->WritePlotFile(plotfilename, "CIC_particles"); } #endif // AMREX_PARTICLES #endif @@ -2297,7 +2505,10 @@ template void AMRSimulation::WriteCheckpointFile // write particle data #ifdef AMREX_PARTICLES if (do_tracers != 0) { - TracerPC->Checkpoint(checkpointname, "particles", true); + TracerPC->Checkpoint(checkpointname, "tracer_particles", true); + } + if (do_cic_particles != 0) { + CICParticles->Checkpoint(checkpointname, "CIC_particles", true); } #endif @@ -2460,7 +2671,12 @@ template void AMRSimulation::ReadCheckpointFile( if (do_tracers != 0) { AMREX_ASSERT(TracerPC == nullptr); TracerPC = std::make_unique(this); - TracerPC->Restart(restart_chkfile, "particles"); + TracerPC->Restart(restart_chkfile, "tracer_particles"); + } + if (do_cic_particles != 0) { + AMREX_ASSERT(CICParticles == nullptr); + CICParticles = std::make_unique(this); + CICParticles->Restart(restart_chkfile, "CIC_particles"); } #endif diff --git a/tests/BinaryOrbit.in b/tests/BinaryOrbit.in new file mode 100644 index 000000000..19bf49341 --- /dev/null +++ b/tests/BinaryOrbit.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 diff --git a/tests/BinaryOrbit_particles.txt b/tests/BinaryOrbit_particles.txt new file mode 100644 index 000000000..e3998ab9d --- /dev/null +++ b/tests/BinaryOrbit_particles.txt @@ -0,0 +1,3 @@ +2 +3.125e12 0. 0. 2.0e34 0. 10332860. 0. +-3.125e12 0. 0. 2.0e34 0. -10332860. 0. diff --git a/tests/SphericalCollapse.in b/tests/SphericalCollapse.in index 943e37065..cb5c328f9 100644 --- a/tests/SphericalCollapse.in +++ b/tests/SphericalCollapse.in @@ -14,16 +14,16 @@ 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.max_level = 1 # 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.max_grid_size = 64 # at least 128 for GPUs amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells -amr.grid_eff = 0.7 # default +amr.grid_eff = 1.0 hydro.reconstruction_order = 3 # PPM cfl = 0.25 max_timesteps = 50000 -stop_time = 0.15 +stop_time = 0.05 derived_vars = gpot @@ -31,7 +31,8 @@ gravity.Gconst = 1.0 # gravitational constant do_reflux = 1 do_subcycle = 0 -do_tracers = 1 # turn on tracer particles +do_tracers = 0 # turn on tracer particles +do_cic_particles = 1 # turns on CIC particles ascent_interval = 100 plotfile_interval = 200