add tracer particles

docs/parameters.rst

@@ -51,7 +51,10 @@ These parameters are read in the `AMRSimulation<problem_t>::readParameters()` fu
      - The time interval (in simulated time) between checkpoint outputs.
    * - do_reflux
      - Integer
-     - this turns on refluxing at coarse-fine boundaries (1) or turns it off (0). Except for debugging, this should always be on when AMR is used.
+     - This turns on refluxing at coarse-fine boundaries (1) or turns it off (0). Except for debugging, this should always be on when AMR is used.
+   * - do_tracers
+     - Integer
+     - This turns on tracer particles. They are initialized one-per-cell and they follow the fluid velocity. Default: 0 (off).
    * - suppress_output
      - Integer
      - If set to 1, this disables output to stdout while the simulation is running.

src/AdvectionSimulation.hpp

@@ -71,6 +71,7 @@ template <typename problem_t> class AdvectionSimulation : public AMRSimulation<p
 	auto computeExtraPhysicsTimestep(int level) -> amrex::Real override;
 	void preCalculateInitialConditions() override;
 	void setInitialConditionsOnGrid(quokka::grid grid_elem) override;
+	void setInitialConditionsOnGridFaceVars(quokka::grid grid_elem) 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;
@@ -147,6 +148,13 @@ template <typename problem_t> void AdvectionSimulation<problem_t>::setInitialCon
 	// user should implement using problem-specific template specialization
 }
 
+template <typename problem_t> void AdvectionSimulation<problem_t>::setInitialConditionsOnGridFaceVars(quokka::grid grid_elem)
+{
+	// default empty implementation
+	// user should implement using problem-specific template specialization
+	// note: an implementation is only required if face-centered vars are used
+}
+
 template <typename problem_t> void AdvectionSimulation<problem_t>::computeAfterTimestep()
 {
 	// do nothing -- user should implement using problem-specific template specialization

src/FCQuantities/test_fc_quantities.cpp

@@ -78,32 +78,33 @@ template <> void RadhydroSimulation<FCQuantities>::setInitialConditionsOnGrid(qu
 	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
 	const amrex::Array4<double> &state = grid_elem.array_;
 	const amrex::Box &indexRange = grid_elem.indexRange_;
-	const quokka::centering cen = grid_elem.cen_;
-	const quokka::direction dir = grid_elem.dir_;
 
-	if (cen == quokka::centering::cc) {
-		const int ncomp_cc = Physics_Indices<FCQuantities>::nvarTotal_cc;
-		// loop over the grid and set the initial condition
-		amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-			for (int n = 0; n < ncomp_cc; ++n) {
-				state(i, j, k, n) = 0; // fill unused quantities with zeros
-			}
-			computeWaveSolution(i, j, k, state, dx, prob_lo);
-		});
-	} else if (cen == quokka::centering::fc) {
-		if (dir == quokka::direction::x) {
-			amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-				state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 1.0 + (i % 2);
-			});
-		} else if (dir == quokka::direction::y) {
-			amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-				state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 2.0 + (j % 2);
-			});
-		} else if (dir == quokka::direction::z) {
-			amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-				state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 3.0 + (k % 2);
-			});
+	const int ncomp_cc = Physics_Indices<FCQuantities>::nvarTotal_cc;
+	// loop over the grid and set the initial condition
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		for (int n = 0; n < ncomp_cc; ++n) {
+			state(i, j, k, n) = 0; // fill unused quantities with zeros
 		}
+		computeWaveSolution(i, j, k, state, dx, prob_lo);
+	});
+}
+
+template <> void RadhydroSimulation<FCQuantities>::setInitialConditionsOnGridFaceVars(quokka::grid grid_elem)
+{
+	// extract grid information
+	const amrex::Array4<double> &state = grid_elem.array_;
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const quokka::direction dir = grid_elem.dir_;
+
+	if (dir == quokka::direction::x) {
+		amrex::ParallelFor(
+		    indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 1.0 + (i % 2); });
+	} else if (dir == quokka::direction::y) {
+		amrex::ParallelFor(
+		    indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 2.0 + (j % 2); });
+	} else if (dir == quokka::direction::z) {
+		amrex::ParallelFor(
+		    indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<FCQuantities>::bfield_index) = 3.0 + (k % 2); });
 	}
 }
 
@@ -159,7 +160,7 @@ auto problem_main() -> int
 	amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state_new_fc_write = sim_write.getNewMF_fc();
 	amrex::Print() << "\n";
 
-	RadhydroSimulation<FCQuantities> sim_restart(BCs_cc);
+	RadhydroSimulation<FCQuantities> sim_restart(BCs_cc, BCs_fc);
 	sim_restart.setChkFile("chk00000");
 	sim_restart.setInitialConditions();
 	amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state_new_fc_restart = sim_restart.getNewMF_fc();

src/RadhydroSimulation.hpp

@@ -13,17 +13,20 @@
 #include <climits>
 #include <filesystem>
 #include <limits>
+#include <memory>
 #include <string>
 #include <tuple>
 #include <unordered_map>
 #include <utility>
 
 #include "AMReX.H"
 #include "AMReX_Algorithm.H"
+#include "AMReX_AmrParticles.H"
 #include "AMReX_Arena.H"
 #include "AMReX_Array.H"
 #include "AMReX_Array4.H"
 #include "AMReX_BCRec.H"
+#include "AMReX_BC_TYPES.H"
 #include "AMReX_BLassert.H"
 #include "AMReX_Box.H"
 #include "AMReX_FArrayBox.H"
@@ -47,6 +50,7 @@
 #include "AMReX_REAL.H"
 #include "AMReX_Utility.H"
 #include "AMReX_YAFluxRegister.H"
+#include "physics_numVars.hpp"
 
 #ifdef AMREX_USE_ASCENT
 #include "AMReX_Conduit_Blueprint.H"
@@ -73,6 +77,9 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	using AMRSimulation<problem_t>::state_old_cc_;
 	using AMRSimulation<problem_t>::state_new_cc_;
 	using AMRSimulation<problem_t>::max_signal_speed_;
+	using AMRSimulation<problem_t>::state_old_fc_;
+	using AMRSimulation<problem_t>::state_new_fc_;
+	using AMRSimulation<problem_t>::TracerPC;
 
 	using AMRSimulation<problem_t>::nghost_cc_;
 	using AMRSimulation<problem_t>::areInitialConditionsDefined_;
@@ -92,6 +99,7 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	using AMRSimulation<problem_t>::finest_level;
 	using AMRSimulation<problem_t>::finestLevel;
 	using AMRSimulation<problem_t>::do_reflux;
+	using AMRSimulation<problem_t>::do_tracers;
 	using AMRSimulation<problem_t>::Verbose;
 	using AMRSimulation<problem_t>::constantDt_;
 	using AMRSimulation<problem_t>::boxArray;
@@ -143,30 +151,19 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	// member functions
 	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc, amrex::Vector<amrex::BCRec> &BCs_fc) : AMRSimulation<problem_t>(BCs_cc, BCs_fc)
 	{
-		defineComponentNames();
-		// read in runtime parameters
-		readParmParse();
-		// set gamma
-		amrex::ParmParse eos("eos");
-		eos.add("eos_gamma", quokka::EOS_Traits<problem_t>::gamma);
-
-		// initialize Microphysics params
-		init_extern_parameters();
-		// initialize Microphysics EOS
-		amrex::Real small_temp = 1e-10;
-		amrex::Real small_dens = 1e-100;
-		eos_init(small_temp, small_dens);
+		initialize();
 	}
 
-	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc) : AMRSimulation<problem_t>(BCs_cc)
+	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc) : AMRSimulation<problem_t>(BCs_cc) { initialize(); }
+
+	inline void initialize()
 	{
 		defineComponentNames();
 		// read in runtime parameters
 		readParmParse();
 		// set gamma
 		amrex::ParmParse eos("eos");
 		eos.add("eos_gamma", quokka::EOS_Traits<problem_t>::gamma);
-
 		// initialize Microphysics params
 		init_extern_parameters();
 		// initialize Microphysics EOS
@@ -184,6 +181,7 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	auto computeExtraPhysicsTimestep(int lev) -> amrex::Real override;
 	void preCalculateInitialConditions() override;
 	void setInitialConditionsOnGrid(quokka::grid grid_elem) override;
+	void setInitialConditionsOnGridFaceVars(quokka::grid grid_elem) override;
 	void advanceSingleTimestepAtLevel(int lev, amrex::Real time, amrex::Real dt_lev, int ncycle) override;
 	void computeAfterTimestep() override;
 	void computeAfterLevelAdvance(int lev, amrex::Real time, amrex::Real dt_lev, int /*ncycle*/);
@@ -313,6 +311,11 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::defineComponen
 	}
 
 	// face-centred
+
+	// add face-centered velocities
+	for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
+		componentNames_fc_.push_back({quokka::face_dir_str[idim] + "-velocity"});
+	}
 	// add mhd state variables
 	if constexpr (Physics_Traits<problem_t>::is_mhd_enabled) {
 		for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
@@ -475,6 +478,13 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::setInitialCond
 	// user should implement using problem-specific template specialization
 }
 
+template <typename problem_t> void RadhydroSimulation<problem_t>::setInitialConditionsOnGridFaceVars(quokka::grid grid_elem)
+{
+	// default empty implementation
+	// user should implement using problem-specific template specialization
+	// note: an implementation is only required if face-centered vars are used
+}
+
 template <typename problem_t> void RadhydroSimulation<problem_t>::computeAfterTimestep()
 {
 	// do nothing -- user should implement if desired
@@ -631,6 +641,9 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::advanceSingleT
 
 	// since we are starting a new timestep, need to swap old and new state vectors
 	std::swap(state_old_cc_[lev], state_new_cc_[lev]);
+	if (Physics_Indices<problem_t>::nvarTotal_fc > 0) {
+		std::swap(state_old_fc_[lev], state_new_fc_[lev]);
+	}
 
 	// check hydro states before update (this can be caused by the flux register!)
 	CHECK_HYDRO_STATES(state_old_cc_[lev]);
@@ -753,9 +766,9 @@ void RadhydroSimulation<problem_t>::FillPatch(int lev, amrex::Real time, amrex::
 	}
 
 	if (cen == quokka::centering::cc) {
-		FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, BCs_cc_, fptype, PreInterpState, PostInterpState);
+		FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, BCs_cc_, cen, fptype, PreInterpState, PostInterpState);
 	} else if (cen == quokka::centering::fc) {
-		FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, BCs_fc_, fptype, PreInterpState, PostInterpState);
+		FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, BCs_fc_, cen, fptype, PreInterpState, PostInterpState);
 	}
 }
 
@@ -857,6 +870,14 @@ void RadhydroSimulation<problem_t>::advanceHydroAtLevelWithRetries(int lev, amre
 		amrex::Copy(originalFineData, fineData, 0, 0, fineData.nComp(), 0);
 	}
 
+#ifdef AMREX_PARTICLES
+	amrex::AmrTracerParticleContainer::ContainerLike<amrex::DefaultAllocator> originalTracerPC;
+	if (do_tracers != 0) {
+		// save the pre-advance tracer particles
+		originalTracerPC = TracerPC->make_alike();
+	}
+#endif
+
 	for (int retry_count = 0; retry_count <= max_retries; ++retry_count) {
 		// reduce timestep by a factor of 2^retry_count
 		const int nsubsteps = std::pow(2, retry_count);
@@ -875,6 +896,13 @@ void RadhydroSimulation<problem_t>::advanceHydroAtLevelWithRetries(int lev, amre
 			if (fr_as_fine != nullptr) {
 				amrex::Copy(fr_as_fine->getFineData(), originalFineData, 0, 0, originalFineData.nComp(), 0);
 			}
+
+#ifdef AMREX_PARTICLES
+			if (do_tracers != 0) {
+				// reset the tracer particles to their pre-advance state
+				TracerPC->copyParticles(originalTracerPC, true);
+			}
+#endif
 		}
 
 		// create temporary multifab for old state
@@ -993,14 +1021,20 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 	amrex::MultiFab state_inter_cc_(grids[lev], dmap[lev], Physics_Indices<problem_t>::nvarTotal_cc, nghost_cc_);
 	state_inter_cc_.setVal(0); // prevent assert in fillBoundaryConditions when radiation is enabled
 
-	// create temporary multifabs for combined RK2 flux
+	// create temporary multifabs for combined RK2 flux and time-average face velocity
 	std::array<amrex::MultiFab, AMREX_SPACEDIM> flux_rk2;
+	std::array<amrex::MultiFab, AMREX_SPACEDIM> avgFaceVel;
+	const int nghost_vel = 2; // 2 ghost faces are needed for tracer particles
 	auto ba = grids[lev];
 	auto dm = dmap[lev];
 	for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
 		auto ba_face = amrex::convert(ba, amrex::IntVect::TheDimensionVector(idim));
+		// initialize flux MultiFab
 		flux_rk2[idim] = amrex::MultiFab(ba_face, dm, ncompHydro_, 0);
 		flux_rk2[idim].setVal(0);
+		// initialize velocity MultiFab
+		avgFaceVel[idim] = amrex::MultiFab(ba_face, dm, 1, nghost_vel);
+		avgFaceVel[idim].setVal(0);
 	}
 
 	// update ghost zones [old timestep]
@@ -1035,6 +1069,7 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 
 		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
 			amrex::MultiFab::Saxpy(flux_rk2[idim], 0.5, fluxArrays[idim], 0, 0, ncompHydro_, 0);
+			amrex::MultiFab::Saxpy(avgFaceVel[idim], 0.5, faceVel[idim], 0, 0, 1, 0);
 		}
 
 		amrex::MultiFab rhs(grids[lev], dmap[lev], ncompHydro_, 0);
@@ -1145,14 +1180,15 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 
 		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
 			amrex::MultiFab::Saxpy(flux_rk2[idim], 0.5, fluxArrays[idim], 0, 0, ncompHydro_, 0);
+			amrex::MultiFab::Saxpy(avgFaceVel[idim], 0.5, faceVel[idim], 0, 0, 1, 0);
 		}
 
 		amrex::MultiFab rhs(grids[lev], dmap[lev], ncompHydro_, 0);
 		amrex::iMultiFab redoFlag(grids[lev], dmap[lev], 1, 1);
 		redoFlag.setVal(quokka::redoFlag::none);
 
 		HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, flux_rk2, dx, ncompHydro_);
-		HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateInter, dx, faceVel, redoFlag);
+		HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateOld, dx, avgFaceVel, redoFlag);
 		HydroSystem<problem_t>::PredictStep(stateOld, stateFinal, rhs, dt_lev, ncompHydro_, redoFlag);
 
 		// do first-order flux correction (FOFC)
@@ -1171,11 +1207,11 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 
 			// replace fluxes around troubled cells with Godunov fluxes
 			replaceFluxes(flux_rk2, FOfluxArrays, redoFlag);
-			replaceFluxes(faceVel, FOfaceVel, redoFlag); // needed for dual energy
+			replaceFluxes(avgFaceVel, FOfaceVel, redoFlag); // needed for dual energy
 
 			// re-do RK update
 			HydroSystem<problem_t>::ComputeRhsFromFluxes(rhs, flux_rk2, dx, ncompHydro_);
-			HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateInter, dx, faceVel, redoFlag);
+			HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateOld, dx, avgFaceVel, redoFlag);
 			HydroSystem<problem_t>::PredictStep(stateOld, stateFinal, rhs, dt_lev, ncompHydro_, redoFlag);
 
 			amrex::Gpu::streamSynchronizeAll(); // just in case
@@ -1213,6 +1249,33 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 	}
 	amrex::Gpu::streamSynchronizeAll();
 
+	// advect tracer particles using avgFaceVel
+#ifdef AMREX_PARTICLES
+	if (do_tracers != 0) {
+		// copy avgFaceVel to state_new_fc_[lev]
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			amrex::Copy(state_new_fc_[lev][idim], avgFaceVel[idim], Physics_Indices<problem_t>::velFirstIndex, 0,
+				    Physics_NumVars::numVelVars_per_dim, nghost_vel);
+		}
+
+		// fill ghost faces
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			fillBoundaryConditions(state_new_fc_[lev][idim], state_new_fc_[lev][idim], lev, time + 0.5 * dt_lev, quokka::centering::fc,
+					       quokka::direction{idim}, AMRSimulation<problem_t>::InterpHookNone, AMRSimulation<problem_t>::InterpHookNone,
+					       FillPatchType::fillpatch_function);
+		}
+
+		// copy state_new_fc_[lev] to avgFaceVel
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			amrex::Copy(avgFaceVel[idim], state_new_fc_[lev][idim], 0, Physics_Indices<problem_t>::velFirstIndex,
+				    Physics_NumVars::numVelVars_per_dim, nghost_vel);
+		}
+
+		// advect particles
+		TracerPC->AdvectWithUmac(avgFaceVel.data(), lev, dt_lev);
+	}
+#endif
+
 	// do Strang split source terms (second half-step)
 	addStrangSplitSourcesWithBuiltin(state_new_cc_[lev], lev, time + dt_lev, 0.5 * dt_lev);
 

src/physics_info.hpp

@@ -2,6 +2,7 @@
 #define PHYSICS_INFO_HPP_
 
 #include "physics_numVars.hpp"
+#include <AMReX.H>
 
 // this struct is specialized by the user application code.
 template <typename problem_t> struct Physics_Traits {
@@ -38,9 +39,11 @@ template <typename problem_t> struct Physics_Indices {
 	static const int pscalarFirstIndex = Physics_NumVars::numHydroVars;
 	static const int radFirstIndex = pscalarFirstIndex + Physics_Traits<problem_t>::numPassiveScalars;
 	// face-centered
-	static const int nvarPerDim_fc = Physics_NumVars::numMHDVars_per_dim * static_cast<int>(Physics_Traits<problem_t>::is_mhd_enabled);
-	static const int nvarTotal_fc = Physics_NumVars::numMHDVars_tot * static_cast<int>(Physics_Traits<problem_t>::is_mhd_enabled);
-	static const int mhdFirstIndex = 0;
+	static const int nvarPerDim_fc = Physics_NumVars::numVelVars_per_dim * static_cast<int>(Physics_Traits<problem_t>::is_hydro_enabled) +
+					 Physics_NumVars::numMHDVars_per_dim * static_cast<int>(Physics_Traits<problem_t>::is_mhd_enabled);
+	static const int nvarTotal_fc = AMREX_SPACEDIM * nvarPerDim_fc;
+	static const int velFirstIndex = 0;
+	static const int mhdFirstIndex = velFirstIndex + Physics_NumVars::numVelVars_per_dim;
 };
 
 #endif // PHYSICS_INFO_HPP_

src/physics_numVars.hpp

@@ -9,7 +9,9 @@ struct Physics_NumVars {
 	static const int numRadVars = 4;
 	// face-centred
 	static const int numMHDVars_per_dim = 1;
+	static const int numVelVars_per_dim = 1;
 	static const int numMHDVars_tot = AMREX_SPACEDIM * numMHDVars_per_dim;
+	static const int numVelVars_tot = AMREX_SPACEDIM * numVelVars_per_dim;
 };
 
 #endif // PHYSICS_NUMVARS_HPP_