add BinaryOrbitCIC test

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,127 @@
+//==============================================================================
+// 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_BC_TYPES.H"
+#include "AMReX_BLassert.H"
+#include "AMReX_Config.H"
+#include "AMReX_FabArrayUtility.H"
+#include "AMReX_MultiFab.H"
+#include "AMReX_ParallelDescriptor.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_Print.H"
+
+#include "RadhydroSimulation.hpp"
+#include "binary_orbit.hpp"
+#include "hydro_system.hpp"
+
+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 <> 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 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(true);
+	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); });
+	}
+}
+
+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 orbital elements
+	// ...
+
+	int 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/CMakeLists.txt

@@ -125,14 +125,10 @@ set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SO
 #endif()
 #link_libraries(QuokkaObj)
 
-add_subdirectory(StarCluster)
-add_subdirectory(PrimordialChem)
-
 add_subdirectory(Advection)
 add_subdirectory(Advection2D)
 add_subdirectory(AdvectionSemiellipse)
 
-
 add_subdirectory(HydroBlast2D)
 add_subdirectory(HydroBlast3D)
 add_subdirectory(HydroContact)
@@ -177,3 +173,6 @@ add_subdirectory(PassiveScalar)
 add_subdirectory(FCQuantities)
 add_subdirectory(SphericalCollapse)
 add_subdirectory(NSCBC)
+add_subdirectory(StarCluster)
+add_subdirectory(PrimordialChem)
+add_subdirectory(BinaryOrbitCIC)

tests/BinaryOrbit.in

@@ -0,0 +1,34 @@
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo     =  -2.5e13 -2.5e13 -2.5e13
+geometry.prob_hi     =   2.5e13  2.5e13  2.5e13
+geometry.is_periodic =  0    0    0
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v              = 1       # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell          = 32 32 32
+amr.max_level       = 0     # number of levels = max_level + 1
+amr.blocking_factor = 32    # grid size must be divisible by this
+amr.max_grid_size   = 32    # at least 128 for GPUs
+amr.n_error_buf     = 3     # minimum 3 cell buffer around tagged cells
+amr.grid_eff        = 1.0
+
+cfl = 0.3
+max_timesteps = 10000
+stop_time = 1.0e7   # s
+
+do_reflux = 1
+do_subcycle = 0
+do_cic_particles = 1    # turns on CIC particles
+
+ascent_interval = -1
+plotfile_interval = 50
+checkpoint_interval = 2000
+derived_vars = gpot

tests/BinaryOrbit_particles.txt

@@ -0,0 +1,3 @@
+2
+3.125e12 0. 0. 2.0e34 0. 10332860. 0.
+-3.125e12 0. 0. 2.0e34 0. -10332860. 0.