add CIC particles (#502)

### Description
This adds cloud-in-cell particles that interact with the gravitational
field and are integrated in time using a kick-drift-kick leapfrog
assuming a global simulation timestep (i.e., no AMR subcycling).

### Related issues
Closes #491.

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

src/AdvectionSimulation.hpp

@@ -72,6 +72,7 @@ template <typename problem_t> class AdvectionSimulation : public AMRSimulation<p
 	void preCalculateInitialConditions() override;
 	void setInitialConditionsOnGrid(quokka::grid grid_elem) override;
 	void setInitialConditionsOnGridFaceVars(quokka::grid grid_elem) override;
+	void createInitialParticles() override;
 	void advanceSingleTimestepAtLevel(int lev, amrex::Real time, amrex::Real dt_lev, int /*ncycle*/) override;
 	void computeAfterTimestep() override;
 	void computeAfterEvolve(amrex::Vector<amrex::Real> &initSumCons) override;
@@ -155,6 +156,13 @@ template <typename problem_t> void AdvectionSimulation<problem_t>::setInitialCon
 	// note: an implementation is only required if face-centered vars are used
 }
 
+template <typename problem_t> void AdvectionSimulation<problem_t>::createInitialParticles()
+{
+	// default empty implementation
+	// user should implement using problem-specific template specialization
+	// note: an implementation is only required if particles are used
+}
+
 template <typename problem_t> void AdvectionSimulation<problem_t>::computeAfterTimestep()
 {
 	// do nothing -- user should implement using problem-specific template specialization

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()

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 <algorithm>
+
+struct BinaryOrbit {
+};
+
+template <> struct quokka::EOS_Traits<BinaryOrbit> {
+	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<BinaryOrbit> {
+	static constexpr bool reconstruct_eint = false;
+};
+
+template <> struct Physics_Traits<BinaryOrbit> {
+	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<BinaryOrbit> {
+	std::vector<amrex::ParticleReal> time{};
+	std::vector<amrex::ParticleReal> dist{};
+};
+
+template <> void RadhydroSimulation<BinaryOrbit>::setInitialConditionsOnGrid(quokka::grid grid_elem)
+{
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const amrex::Array4<double> &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<BinaryOrbit>::density_index) = rho;
+		state_cc(i, j, k, HydroSystem<BinaryOrbit>::x1Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<BinaryOrbit>::x2Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<BinaryOrbit>::x3Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<BinaryOrbit>::energy_index) = 0;
+		state_cc(i, j, k, HydroSystem<BinaryOrbit>::internalEnergy_index) = 0;
+	});
+}
+
+template <> void RadhydroSimulation<BinaryOrbit>::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<BinaryOrbit>::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<BinaryOrbit>::computeAfterTimestep()
+{
+	// every N cycles, save particle statistics
+	static int cycle = 1;
+	if (cycle % 10 == 0) {
+		// create single-box particle container
+		amrex::ParticleContainer<quokka::CICParticleRealComps> 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<quokka::CICParticleContainer::ParticleType> 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<BinaryOrbit>::x1Momentum_index) && (dim == 0)) {
+			return true;
+		}
+		if ((n == HydroSystem<BinaryOrbit>::x2Momentum_index) && (dim == 1)) {
+			return true;
+		}
+		if ((n == HydroSystem<BinaryOrbit>::x3Momentum_index) && (dim == 2)) {
+			return true;
+		}
+		return false;
+	};
+
+	const int ncomp_cc = Physics_Indices<BinaryOrbit>::nvarTotal_cc;
+	amrex::Vector<amrex::BCRec> 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<BinaryOrbit> 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<double>::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;
+}

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 <fstream>
+
+// internal headers
+#include "hydro_system.hpp"
+#include "interpolate.hpp"
+
+// function definitions
+auto problem_main() -> int;
+
+#endif // TEST_BINARY_ORBIT_HPP_

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)

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<CICParticleRealComps>;
+using CICParticleIterator = amrex::ParIter<CICParticleRealComps>;
+
+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<amrex::Real> const &rho,
+							    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &plo,
+							    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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_

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)