add CIC deposition

src/CICParticles.hpp

@@ -9,17 +9,40 @@
 /// \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_ParticleInterpolators.H"
+#include "AMReX_Particles.H"
 namespace quokka
 {
 
 constexpr int CICParticleReals = 4; // mass vx vy vz
 constexpr int CICParticleInts = 0;
 
-enum ParticleDataIdx {ParticleMass = 0, ParticleVx, ParticleVy, ParticleVz};
+enum ParticleDataIdx { ParticleMassIdx = 0, ParticleVxIdx, ParticleVyIdx, ParticleVzIdx };
 
 using CICParticleContainer = amrex::AmrParticleContainer<0, 0, CICParticleReals, CICParticleInts>;
+
+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/RadhydroSimulation.hpp

@@ -110,6 +110,7 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	using AMRSimulation<problem_t>::WriteCheckpointFile;
 	using AMRSimulation<problem_t>::GetData;
 	using AMRSimulation<problem_t>::FillPatchWithData;
+	using AMRSimulation<problem_t>::Gconst_;
 
 	using AMRSimulation<problem_t>::densityFloor_;
 	using AMRSimulation<problem_t>::tempFloor_;
@@ -146,8 +147,6 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 
 	amrex::Long radiationCellUpdates_ = 0; // total number of radiation cell-updates
 
-	amrex::Real Gconst_ = C::Gconst; // gravitational constant G
-
 	// member functions
 	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc, amrex::Vector<amrex::BCRec> &BCs_fc) : AMRSimulation<problem_t>(BCs_cc, BCs_fc)
 	{
@@ -353,12 +352,6 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::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");

src/simulation.hpp

@@ -72,11 +72,12 @@ namespace filesystem = experimental::filesystem;
 #include <AMReX_Utility.H>
 #include <fmt/core.h>
 #include <yaml-cpp/yaml.h>
+#include "fundamental_constants.H"
 
 #ifdef AMREX_PARTICLES
+#include "CICParticles.hpp"
 #include <AMReX_AmrParticles.H>
 #include <AMReX_Particles.H>
-#include "CICParticles.hpp"
 #endif
 
 #if AMREX_SPACEDIM == 3
@@ -386,9 +387,12 @@ template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 	amrex::Long cellUpdates_ = 0;
 	amrex::Vector<amrex::Long> 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;
@@ -622,6 +626,12 @@ template <typename problem_t> void AMRSimulation<problem_t>::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
+	{
+		amrex::ParmParse hpp("gravity");
+		hpp.query("Gconst", Gconst_);
+	}
 }
 
 template <typename problem_t> void AMRSimulation<problem_t>::setInitialConditions()
@@ -1004,12 +1014,14 @@ template <typename problem_t> void AMRSimulation<problem_t>::calculateGpotAllLev
 			rhs[lev].define(grids[lev], dmap[lev], ncomp, nghost);
 			phi[lev].setVal(0); // set initial guess to zero
 			rhs[lev].setVal(0);
-			fillPoissonRhsAtLevel(rhs[lev], lev);
 		}
 
-		// TODO(bwibking): deposit particles using amrex::ParticleToMesh
+		// deposit particles using amrex::ParticleToMesh
+		amrex::ParticleToMesh(*CICParticles, amrex::GetVecOfPtrs(rhs), 0, finest_level, quokka::CICDeposition{Gconst_, quokka::ParticleMassIdx, 0, 1});
 
+		// add fluid density
 		for (int lev = 0; lev <= finest_level; ++lev) {
+			fillPoissonRhsAtLevel(rhs[lev], lev);
 			rhs_min = std::min(rhs_min, rhs[lev].min(0));
 		}