updated

.github/workflows/clang-tidy-comments.yml

@@ -14,7 +14,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: ZedThree/clang-tidy-review/post@v0.15.1
+      - uses: ZedThree/clang-tidy-review/post@v0.17.0
       # lgtm_comment_body, max_comments, and annotations need to be set on the posting workflow in a split setup
         with:
         # adjust options as necessary

.github/workflows/clang-tidy.yml

@@ -17,7 +17,7 @@ jobs:
         submodules: true
         fetch-depth: 0
 
-    - uses: ZedThree/clang-tidy-review@v0.15.1
+    - uses: ZedThree/clang-tidy-review@v0.17.0
       id: review
       with:
         config_file: src/.clang-tidy
@@ -27,7 +27,7 @@ jobs:
         split_workflow: true
 
     # Uploads an artefact containing clang_fixes.json
-    - uses: ZedThree/clang-tidy-review/upload@v0.15.1
+    - uses: ZedThree/clang-tidy-review/upload@v0.17.0
 
     # If there are any comments, fail the check
     - if: steps.review.outputs.total_comments > 0

.github/workflows/codespell.yml

@@ -24,7 +24,7 @@ jobs:
           python-version: '3.10'
 
       - name: Cache pip
-        uses: actions/cache@v3
+        uses: actions/cache@v4
         with:
           # this path is specific to Ubuntu
           path: ~/.cache/pip

.github/workflows/hip.yml

@@ -28,7 +28,7 @@ jobs:
     - name: Python and HDF5 dependencies
       run: sudo apt-get update && sudo apt-get install python3-dev python3-numpy python3-matplotlib python3-pip libhdf5-mpi-dev
     - name: Set Up Cache
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       with:
         path: ~/.cache/ccache
         key: ccache-${{ github.workflow }}-${{ github.job }}-git-${{ github.sha }}

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

@@ -179,3 +179,6 @@ add_subdirectory(PassiveScalar)
 add_subdirectory(FCQuantities)
 add_subdirectory(SphericalCollapse)
 add_subdirectory(NSCBC)
+add_subdirectory(StarCluster)
+add_subdirectory(PrimordialChem)
+add_subdirectory(BinaryOrbitCIC)