[WIP] add periodic Poisson gravity BC

src/problems/Jeans/CMakeLists.txt

@@ -0,0 +1,8 @@
+if (AMReX_SPACEDIM EQUAL 3)
+    add_executable(jeans_collapse jeans.cpp ${QuokkaObjSources})
+    if(AMReX_GPU_BACKEND MATCHES "CUDA")
+        setup_target_for_cuda_compilation(jeans_collapse)
+    endif()
+
+    add_test(NAME JeansCollapse COMMAND jeans_collapse Jeans.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()

src/problems/Jeans/jeans.cpp

@@ -0,0 +1,166 @@
+//==============================================================================
+// 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 jeans.cpp
+/// \brief Defines a test problem for Jeans collapse.
+///
+
+#include "AMReX_BC_TYPES.H"
+#include "AMReX_BLassert.H"
+#include "AMReX_MultiFab.H"
+#include "AMReX_ParmParse.H"
+
+#include "RadhydroSimulation.hpp"
+#include "hydro/hydro_system.hpp"
+#include "jeans.hpp"
+
+struct CollapseProblem {
+};
+
+template <> struct quokka::EOS_Traits<CollapseProblem> {
+	static constexpr double gamma = 5. / 3.;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
+};
+
+template <> struct HydroSystem_Traits<CollapseProblem> {
+	static constexpr bool reconstruct_eint = false;
+};
+
+template <> struct Physics_Traits<CollapseProblem> {
+	// cell-centred
+	static constexpr bool is_hydro_enabled = true;
+	static constexpr int numMassScalars = 0;		     // number of mass scalars
+	static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars
+	static constexpr bool is_radiation_enabled = false;
+	// face-centred
+	static constexpr bool is_mhd_enabled = false;
+	static constexpr int nGroups = 1; // number of radiation groups
+};
+
+template <> void RadhydroSimulation<CollapseProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)
+{
+	// set initial conditions
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = grid_elem.dx_;
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_hi = grid_elem.prob_hi_;
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const amrex::Array4<double> &state_cc = grid_elem.array_;
+
+	amrex::Real x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]);
+	amrex::Real y0 = prob_lo[1] + 0.5 * (prob_hi[1] - prob_lo[1]);
+	amrex::Real z0 = prob_lo[2] + 0.5 * (prob_hi[2] - prob_lo[2]);
+
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		amrex::Real const x = prob_lo[0] + (i + static_cast<amrex::Real>(0.5)) * dx[0];
+		amrex::Real const y = prob_lo[1] + (j + static_cast<amrex::Real>(0.5)) * dx[1];
+		amrex::Real const z = prob_lo[2] + (k + static_cast<amrex::Real>(0.5)) * dx[2];
+		amrex::Real const r = std::sqrt(std::pow(x - x0, 2) + std::pow(y - y0, 2) + std::pow(z - z0, 2));
+
+		double rho_min = 1.0e-5;
+		double rho_max = 10.0;
+		double R_sphere = 0.5;
+		double R_smooth = 0.025;
+		double rho = std::max(rho_min, rho_max * ((std::tanh((R_sphere - r) / R_smooth) + 1.0) / 2.0));
+		double P = 1.0e-1;
+
+		AMREX_ASSERT(!std::isnan(rho));
+		AMREX_ASSERT(!std::isnan(P));
+
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::density_index) = rho;
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::x1Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::x2Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::x3Momentum_index) = 0;
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::energy_index) = quokka::EOS<CollapseProblem>::ComputeEintFromPres(rho, P);
+		state_cc(i, j, k, HydroSystem<CollapseProblem>::internalEnergy_index) = quokka::EOS<CollapseProblem>::ComputeEintFromPres(rho, P);
+	});
+}
+
+template <> void RadhydroSimulation<CollapseProblem>::createInitialParticles()
+{
+	// add particles at random positions in the box
+	const bool generate_on_root_rank = true;
+	const int iseed = 42;
+	const int num_particles = 1000;
+	const double total_particle_mass = 0.5; // about 0.1 of the total fluid mass
+	const double particle_mass = total_particle_mass / static_cast<double>(num_particles);
+
+	const quokka::CICParticleContainer::ParticleInitData pdata = {{particle_mass, 0, 0, 0}, {}, {}, {}}; // {mass vx vy vz}, empty, empty, empty
+	CICParticles->InitRandom(num_particles, iseed, pdata, generate_on_root_rank);
+}
+
+template <> void RadhydroSimulation<CollapseProblem>::ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real /*time*/, int /*ngrow*/)
+{
+	// tag cells for refinement
+	const Real q_min = 5.0; // minimum density for refinement
+
+	for (amrex::MFIter mfi(state_new_cc_[lev]); mfi.isValid(); ++mfi) {
+		const amrex::Box &box = mfi.validbox();
+		const auto state = state_new_cc_[lev].const_array(mfi);
+		const auto tag = tags.array(mfi);
+		const int nidx = HydroSystem<CollapseProblem>::density_index;
+
+		amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+			Real const q = state(i, j, k, nidx);
+			if (q > q_min) {
+				tag(i, j, k) = amrex::TagBox::SET;
+			}
+		});
+	}
+}
+
+template <> void RadhydroSimulation<CollapseProblem>::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); });
+	}
+}
+
+auto problem_main() -> int
+{
+	auto isNormalComp = [=](int n, int dim) {
+		if ((n == HydroSystem<CollapseProblem>::x1Momentum_index) && (dim == 0)) {
+			return true;
+		}
+		if ((n == HydroSystem<CollapseProblem>::x2Momentum_index) && (dim == 1)) {
+			return true;
+		}
+		if ((n == HydroSystem<CollapseProblem>::x3Momentum_index) && (dim == 2)) {
+			return true;
+		}
+		return false;
+	};
+
+	const int ncomp_cc = Physics_Indices<CollapseProblem>::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<CollapseProblem> sim(BCs_cc);
+	sim.doPoissonSolve_ = 1; // enable self-gravity
+
+	// initialize
+	sim.setInitialConditions();
+
+	// evolve
+	sim.evolve();
+
+	int status = 0;
+	return status;
+}

src/problems/Jeans/jeans.hpp

@@ -0,0 +1,22 @@
+#ifndef TEST_SPHERICAL_COLLAPSE_HPP_ // NOLINT
+#define TEST_SPHERICAL_COLLAPSE_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 spherical_collapse.hpp
+/// \brief Defines a test problem for a pressureless spherical collapse
+///
+
+// external headers
+#include <fstream>
+
+// internal headers
+#include "hydro/hydro_system.hpp"
+#include "math/interpolate.hpp"
+
+// function definitions
+auto problem_main() -> int;
+
+#endif // TEST_SPHERICAL_COLLAPSE_HPP_

src/simulation.hpp

@@ -98,8 +98,6 @@ using namespace conduit;
 using namespace ascent;
 #endif
 
-enum class ParticleStep { BeforePoissonSolve, AfterPoissonSolve };
-
 using variant_t = std::variant<amrex::Real, std::string>;
 
 namespace YAML
@@ -137,6 +135,8 @@ template <> struct as_if<std::string, std::optional<std::string>> {
 
 enum class FillPatchType { fillpatch_class, fillpatch_function };
 
+enum class GravityBCType { periodic, open };
+
 // Main simulation class; solvers should inherit from this
 template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 {
@@ -154,19 +154,20 @@ template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 	amrex::Real cflNumber_ = 0.3;	      // default
 	amrex::Real dtToleranceFactor_ = 1.1; // default
 	amrex::Long cycleCount_ = 0;
-	amrex::Long maxTimesteps_ = 1e4;	    // default
-	amrex::Long maxWalltime_ = 0;		    // default: no limit
-	int ascentInterval_ = -1;		    // -1 == no in-situ renders with Ascent
-	int plotfileInterval_ = -1;		    // -1 == no output
-	int projectionInterval_ = -1;		    // -1 == no output
-	int statisticsInterval_ = -1;		    // -1 == no output
-	amrex::Real plotTimeInterval_ = -1.0;	    // time interval for plt file
-	amrex::Real checkpointTimeInterval_ = -1.0; // time interval for checkpoints
-	int checkpointInterval_ = -1;		    // -1 == no output
-	int amrInterpMethod_ = 1;		    // 0 == piecewise constant, 1 == lincc_interp
-	amrex::Real reltolPoisson_ = 1.0e-5;	    // default
-	amrex::Real abstolPoisson_ = 1.0e-5;	    // default (scaled by minimum RHS value)
-	int doPoissonSolve_ = 0;		    // 1 == self-gravity enabled, 0 == disabled
+	amrex::Long maxTimesteps_ = 1e4;		// default
+	amrex::Long maxWalltime_ = 0;			// default: no limit
+	int ascentInterval_ = -1;			// -1 == no in-situ renders with Ascent
+	int plotfileInterval_ = -1;			// -1 == no output
+	int projectionInterval_ = -1;			// -1 == no output
+	int statisticsInterval_ = -1;			// -1 == no output
+	amrex::Real plotTimeInterval_ = -1.0;		// time interval for plt file
+	amrex::Real checkpointTimeInterval_ = -1.0;	// time interval for checkpoints
+	int checkpointInterval_ = -1;			// -1 == no output
+	int amrInterpMethod_ = 1;			// 0 == piecewise constant, 1 == lincc_interp
+	amrex::Real reltolPoisson_ = 1.0e-5;		// default
+	amrex::Real abstolPoisson_ = 1.0e-5;		// default (scaled by minimum RHS value)
+	GravityBCType gravityBC_ = GravityBCType::open; // default BC
+	int doPoissonSolve_ = 0;			// 1 == self-gravity enabled, 0 == disabled
 	amrex::Vector<amrex::MultiFab> phi;
 
 	amrex::Real densityFloor_ = 0.0; // default
@@ -490,9 +491,8 @@ template <typename problem_t> void AMRSimulation<problem_t>::PerformanceHints()
 		const amrex::Long nboxes = boxArray(ilev).size();
 		if (amrex::ParallelDescriptor::NProcs() > nboxes) {
 			amrex::Print() << "\n[Warning] [Performance] Too many resources / too little work!\n"
-				       << "  It looks like you requested more compute resources than "
-				       << "  the number of boxes of cells available on level " << ilev << " (" << nboxes << "). "
-				       << "You started with (" << amrex::ParallelDescriptor::NProcs() << ") MPI ranks, so ("
+				       << "  It looks like you requested more compute resources than " << "  the number of boxes of cells available on level "
+				       << ilev << " (" << nboxes << "). " << "You started with (" << amrex::ParallelDescriptor::NProcs() << ") MPI ranks, so ("
 				       << amrex::ParallelDescriptor::NProcs() - nboxes << ") rank(s) will have no work on this level.\n"
 #ifdef AMREX_USE_GPU
 				       << "  On GPUs, consider using 1-8 boxes per GPU per level that "
@@ -1018,16 +1018,34 @@ template <typename problem_t> void AMRSimulation<problem_t>::calculateGpotAllLev
 
 		BL_PROFILE_REGION("GravitySolver");
 
-		// set up elliptic solve object
-		amrex::OpenBCSolver poissonSolver(Geom(0, finest_level), boxArray(0, finest_level), DistributionMap(0, finest_level));
-		if (verbose) {
-			poissonSolver.setVerbose(true);
-			poissonSolver.setBottomVerbose(false);
-			amrex::Print() << "Doing Poisson solve...\n\n";
+		// set up Poisson solver
+		std::unique_ptr<amrex::OpenBCSolver> poissonSolverOpenBC;
+		std::unique_ptr<amrex::MLPoisson> poissonSolverPeriodic;
+
+		if (gravityBC_ == GravityBCType::periodic) {
+			// periodic boundary conditions
+
+			poissonSolverPeriodic =
+			    std::make_unique<amrex::MLPoisson>(Geom(0, finest_level), boxArray(0, finest_level), DistributionMap(0, finest_level));
+
+			amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> bc_lo{amrex::LinOpBCType::Periodic};
+			amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> bc_hi{amrex::LinOpBCType::Periodic};
+			poissonSolverPeriodic->setDomainBC(bc_lo, bc_hi);
+		} else {
+			// open boundary conditions
+			// solve Poisson equation with open b.c. using the method of James (1977)
+
+			poissonSolverOpenBC =
+			    std::make_unique<amrex::OpenBCSolver>(Geom(0, finest_level), boxArray(0, finest_level), DistributionMap(0, finest_level));
+
+			if (verbose) {
+				poissonSolverOpenBC->setVerbose(true);
+				poissonSolverOpenBC->setBottomVerbose(false);
+				amrex::Print() << "Doing Poisson solve...\n\n";
+			}
 		}
 
 		phi.resize(finest_level + 1);
-		// solve Poisson equation with open b.c. using the method of James (1977)
 		amrex::Vector<amrex::MultiFab> rhs(finest_level + 1);
 		const int nghost = 1;
 		const int ncomp = 1;
@@ -1055,7 +1073,17 @@ template <typename problem_t> void AMRSimulation<problem_t>::calculateGpotAllLev
 		}
 
 		amrex::Real abstol = abstolPoisson_ * rhs_min;
-		poissonSolver.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), reltolPoisson_, abstol);
+
+		if (gravityBC_ == GravityBCType::periodic) {
+			amrex::MLMG mlmg(*poissonSolverPeriodic);
+			if (verbose) {
+				mlmg.setVerbose(true);
+			}
+			mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), reltolPoisson_, abstol);
+		} else {
+			poissonSolverOpenBC->solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), reltolPoisson_, abstol);
+		}
+
 		if (verbose) {
 			amrex::Print() << "\n";
 		}
@@ -1087,9 +1115,7 @@ template <typename problem_t> void AMRSimulation<problem_t>::ellipticSolveAllLev
 {
 #if AMREX_SPACEDIM == 3
 	if (doPoissonSolve_) {
-
 		calculateGpotAllLevels();
-
 		gravAccelAllLevels(dt);
 	}
 #endif