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