add tracer particles

src/simulation.hpp

@@ -336,6 +336,13 @@ template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 	amrex::Long cellUpdates_ = 0;
 	amrex::Vector<amrex::Long> cellUpdatesEachLevel_;
 
+	// tracer particles
+#ifdef AMREX_PARTICLES
+	void InitParticles(); // create tracer particles
+	int do_tracers = 0;
+	std::unique_ptr<amrex::AmrTracerParticleContainer> TracerPC;
+#endif
+
 	// external objects
 #ifdef AMREX_USE_ASCENT
 	Ascent ascent_;
@@ -568,6 +575,12 @@ template <typename problem_t> void AMRSimulation<problem_t>::setInitialCondition
 		InitFromScratch(time);
 		AverageDown();
 
+#ifdef AMREX_PARTICLES
+		if (do_tracers) {
+			InitParticles();
+		}
+#endif
+
 		if (checkpointInterval_ > 0) {
 			WriteCheckpointFile();
 		}
@@ -976,6 +989,13 @@ template <typename problem_t> auto AMRSimulation<problem_t>::timeStepWithSubcycl
 					}
 				}
 
+#ifdef AMREX_PARTICLES
+				// redistribute particles
+				if (do_tracers) {
+					TracerPC->Redistribute(lev);
+				}
+#endif
+
 				// do fix-up on all levels that have been re-gridded
 				for (int k = lev; k <= finest_level; ++k) {
 					FixupState(k);
@@ -1057,6 +1077,8 @@ template <typename problem_t> auto AMRSimulation<problem_t>::timeStepWithSubcycl
 	// do hyperbolic advance over all levels
 	advanceSingleTimestepAtLevel(lev, time, dt_[lev], nsubsteps[lev]);
 
+	// TODO(bwibking): advect tracer particles here
+
 	++istep[lev];
 	cellUpdates_ += CountCells(lev); // keep track of total number of cell updates
 	cellUpdatesEachLevel_[lev] += CountCells(lev);
@@ -1097,6 +1119,21 @@ template <typename problem_t> auto AMRSimulation<problem_t>::timeStepWithSubcycl
 		fillpatcher_[lev + 1].reset(); // because the data on lev have changed.
 	}
 
+#ifdef AMREX_PARTICLES
+	// redistribute particles
+	if (do_tracers) {
+		int redistribute_ngrow = 0;
+		if ((iteration < nsubsteps[lev]) || (lev == 0)) {
+			if (lev == 0) {
+				redistribute_ngrow = 0;
+			} else {
+				redistribute_ngrow = iteration;
+			}
+			TracerPC->Redistribute(lev, TracerPC->finestLevel(), redistribute_ngrow);
+		}
+	}
+#endif
+
 	return stepsLeft;
 }
 
@@ -1640,6 +1677,22 @@ template <typename problem_t> void AMRSimulation<problem_t>::AverageDownTo(int c
 	}
 }
 
+#ifdef AMREX_PARTICLES
+template <typename problem_t> void AMRSimulation<problem_t>::InitParticles()
+{
+	if (do_tracers) {
+		AMREX_ASSERT(TracerPC == nullptr);
+		TracerPC = std::make_unique<AmrTracerParticleContainer>(this);
+
+		AmrTracerParticleContainer::ParticleInitData pdata = {{AMREX_D_DECL(0.0, 0.0, 0.0)}, {}, {}, {}};
+
+		TracerPC->SetVerbose(0);
+		TracerPC->InitOnePerCell(0.5, 0.5, 0.5, pdata);
+		TracerPC->Redistribute();
+	}
+}
+#endif
+
 // get plotfile name
 template <typename problem_t> auto AMRSimulation<problem_t>::PlotFileName(int lev) const -> std::string { return amrex::Concatenate(plot_file, lev, 5); }
 
@@ -1852,11 +1905,18 @@ template <typename problem_t> void AMRSimulation<problem_t>::WritePlotFile()
 	amrex::Print() << "Writing plotfile " << plotfilename << "\n";
 
 #ifdef QUOKKA_USE_OPENPMD
+	// TODO(bwibking): write particles using openPMD
 	quokka::OpenPMDOutput::WriteFile(varnames, finest_level + 1, mf_ptr, Geom(), plot_file, tNew_[0], istep[0]);
 	WriteMetadataFile(plotfilename + ".yaml");
 #else
 	amrex::WriteMultiLevelPlotfile(plotfilename, finest_level + 1, mf_ptr, varnames, Geom(), tNew_[0], istep, refRatio());
 	WriteMetadataFile(plotfilename + "/metadata.yaml");
+#ifdef AMREX_PARTICLES
+	// write particles
+	if (do_tracers) {
+		TracerPC->WritePlotFile(plotfilename, "particles");
+	}
+#endif // AMREX_PARTICLES
 #endif
 }
 
@@ -2178,6 +2238,13 @@ template <typename problem_t> void AMRSimulation<problem_t>::WriteCheckpointFile
 		}
 	}
 
+	// write particle data
+#ifdef AMREX_PARTICLES
+	if (do_tracers) {
+		TracerPC->Checkpoint(checkpointname, "particles", true);
+	}
+#endif
+
 	// create symlink and point it at this checkpoint dir
 	SetLastCheckpointSymlink(checkpointname);
 }
@@ -2331,6 +2398,16 @@ template <typename problem_t> void AMRSimulation<problem_t>::ReadCheckpointFile(
 			}
 		}
 	}
+
+#ifdef AMREX_PARTICLES
+	// read particle data
+	if (do_tracers) {
+		AMREX_ASSERT(TracerPC == nullptr);
+		TracerPC = std::make_unique<AmrTracerParticleContainer>(this);
+		TracerPC->Restart(restart_chkfile, "particles");
+	}
+#endif
+
 	areInitialConditionsDefined_ = true;
 }