add Navier-Stokes Characteristic Boundary Conditions

src/RadhydroSimulation.hpp

@@ -10,7 +10,10 @@
 /// timestepping, solving, and I/O of a simulation for radiation moments.
 
 #include <array>
+#include <ascent.hpp>
 #include <climits>
+#include <conduit_node.hpp>
+#include <filesystem>
 #include <limits>
 #include <string>
 #include <tuple>
@@ -24,6 +27,7 @@
 #include "AMReX_BCRec.H"
 #include "AMReX_BLassert.H"
 #include "AMReX_Box.H"
+#include "AMReX_Conduit_Blueprint.H"
 #include "AMReX_FArrayBox.H"
 #include "AMReX_FabArray.H"
 #include "AMReX_FabArrayUtility.H"
@@ -36,8 +40,10 @@
 #include "AMReX_IntVect.H"
 #include "AMReX_MultiFab.H"
 #include "AMReX_MultiFabUtil.H"
+#include "AMReX_ParallelDescriptor.H"
 #include "AMReX_ParmParse.H"
 #include "AMReX_PhysBCFunct.H"
+#include "AMReX_PlotFileUtil.H"
 #include "AMReX_Print.H"
 #include "AMReX_REAL.H"
 #include "AMReX_Utility.H"
@@ -70,9 +76,11 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	using AMRSimulation<problem_t>::componentNames_fc_;
 	using AMRSimulation<problem_t>::cflNumber_;
 	using AMRSimulation<problem_t>::fillBoundaryConditions;
+	using AMRSimulation<problem_t>::CustomPlotFileName;
 	using AMRSimulation<problem_t>::geom;
 	using AMRSimulation<problem_t>::grids;
 	using AMRSimulation<problem_t>::dmap;
+	using AMRSimulation<problem_t>::istep;
 	using AMRSimulation<problem_t>::flux_reg_;
 	using AMRSimulation<problem_t>::incrementFluxRegisters;
 	using AMRSimulation<problem_t>::finest_level;
@@ -114,6 +122,7 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	amrex::Real errorNorm_ = NAN;
 	amrex::Real pressureFloor_ = 0.;
 
+	int lowLevelDebuggingOutput_ = 0;	// 0 == do nothing; 1 == output intermediate multifabs used in hydro each timestep (ONLY USE FOR DEBUGGING)
 	int integratorOrder_ = 2;		// 1 == forward Euler; 2 == RK2-SSP (default)
 	int reconstructionOrder_ = 3;		// 1 == donor cell; 2 == PLM; 3 == PPM (default)
 	int radiationReconstructionOrder_ = 3;	// 1 == donor cell; 2 == PLM; 3 == PPM (default)
@@ -310,6 +319,8 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::readParmParse(
 	// set hydro runtime parameters
 	{
 		amrex::ParmParse hpp("hydro");
+		hpp.query("low_level_debugging_output", lowLevelDebuggingOutput_);
+		hpp.query("rk_integrator_order", integratorOrder_);
 		hpp.query("reconstruction_order", reconstructionOrder_);
 		hpp.query("use_dual_energy", useDualEnergy_);
 		hpp.query("abort_on_fofc_failure", abortOnFofcFailure_);
@@ -325,11 +336,14 @@ template <typename problem_t> void RadhydroSimulation<problem_t>::readParmParse(
 	// set cooling runtime parameters
 	{
 		amrex::ParmParse hpp("cooling");
+		int alwaysReadTables = 0;
 		hpp.query("enabled", enableCooling_);
+		hpp.query("read_tables_even_if_disabled", alwaysReadTables);
 		hpp.query("grackle_data_file", coolingTableFilename_);
 
-		if (enableCooling_ == 1) {
+		if ((enableCooling_ == 1) || (alwaysReadTables == 1)) {
 			// read Cloudy tables
+			amrex::Print() << "Reading Cloudy tables...\n";
 			quokka::cooling::readCloudyData(coolingTableFilename_, cloudyTables_);
 		}
 	}
@@ -901,6 +915,18 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 		fillBoundaryConditions(state_old_cc_tmp, state_old_cc_tmp, lev, time, quokka::centering::cc, quokka::direction::na, PreInterpState,
 				       PostInterpState);
 
+		// LOW LEVEL DEBUGGING: output state_old_cc_tmp (with ghost cells)
+		if (lowLevelDebuggingOutput_ == 1) {
+			// write Blueprint HDF5 files
+			conduit::Node mesh;
+			amrex::SingleLevelToBlueprint(state_old_cc_tmp, componentNames_cc_, geom[lev], time, istep[lev] + 1, mesh);
+			amrex::WriteBlueprintFiles(mesh, "debug_stage1_filled_state_old", istep[lev] + 1, "hdf5");
+
+			// write AMReX plotfile
+			// WriteSingleLevelPlotfile(CustomPlotFileName("debug_stage1_filled_state_old", istep[lev]+1),
+			//	state_old_cc_tmp, componentNames_cc_, geom[lev], time, istep[lev]+1);
+		}
+
 		// check state validity
 		AMREX_ASSERT(!state_old_cc_tmp.contains_nan(0, state_old_cc_tmp.nComp()));
 		AMREX_ASSERT(!state_old_cc_tmp.contains_nan()); // check ghost cells
@@ -918,6 +944,32 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 		HydroSystem<problem_t>::AddInternalEnergyPdV(rhs, stateOld, dx, faceVel, redoFlag);
 		HydroSystem<problem_t>::PredictStep(stateOld, stateNew, rhs, dt_lev, ncompHydro_, redoFlag);
 
+		// LOW LEVEL DEBUGGING: output rhs
+		if (lowLevelDebuggingOutput_ == 1) {
+			// write rhs
+			std::string plotfile_name = CustomPlotFileName("debug_stage1_rhs_fluxes", istep[lev] + 1);
+			WriteSingleLevelPlotfile(plotfile_name, rhs, componentNames_cc_, geom[lev], time, istep[lev] + 1);
+
+			// write fluxes
+			for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+				if (amrex::ParallelDescriptor::IOProcessor()) {
+					std::filesystem::create_directories(plotfile_name + "/raw_fields/Level_" + std::to_string(lev));
+				}
+				std::string fullprefix =
+				    amrex::MultiFabFileFullPrefix(lev, plotfile_name, "raw_fields/Level_", std::string("Flux_") + quokka::face_dir_str[idim]);
+				amrex::VisMF::Write(fluxArrays[idim], fullprefix);
+			}
+			// write face velocities
+			for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+				if (amrex::ParallelDescriptor::IOProcessor()) {
+					std::filesystem::create_directories(plotfile_name + "/raw_fields/Level_" + std::to_string(lev));
+				}
+				std::string fullprefix = amrex::MultiFabFileFullPrefix(lev, plotfile_name, "raw_fields/Level_",
+										       std::string("FaceVel_") + quokka::face_dir_str[idim]);
+				amrex::VisMF::Write(faceVel[idim], fullprefix);
+			}
+		}
+
 		// do first-order flux correction (FOFC)
 		amrex::Gpu::streamSynchronizeAll(); // just in case
 		int ncells_bad = redoFlag.sum(0);
@@ -965,7 +1017,7 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 	amrex::Gpu::streamSynchronizeAll();
 
 	// Stage 2 of RK2-SSP
-	{
+	if (integratorOrder_ == 2) {
 		// update ghost zones [intermediate stage stored in state_inter_cc_]
 		fillBoundaryConditions(state_inter_cc_, state_inter_cc_, lev, time + dt_lev, quokka::centering::cc, quokka::direction::na, PreInterpState,
 				       PostInterpState);
@@ -1030,6 +1082,8 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 			// increment flux registers
 			incrementFluxRegisters(fr_as_crse, fr_as_fine, fluxArrays, lev, fluxScaleFactor * dt_lev);
 		}
+	} else { // we are only doing forward Euler
+		amrex::Copy(state_new_cc_[lev], state_inter_cc_, 0, 0, ncompHydro_, 0);
 	}
 	amrex::Gpu::streamSynchronizeAll();
 
@@ -1121,7 +1175,7 @@ auto RadhydroSimulation<problem_t>::computeHydroFluxes(amrex::MultiFab const &co
 	auto ba = grids[lev];
 	auto dm = dmap[lev];
 	const int flatteningGhost = 2;
-	const int reconstructRange = 1;
+	const int reconstructGhost = 1;
 
 	// allocate temporary MultiFabs
 	amrex::MultiFab primVar(ba, dm, nvars, nghost_cc_);
@@ -1137,8 +1191,8 @@ auto RadhydroSimulation<problem_t>::computeHydroFluxes(amrex::MultiFab const &co
 
 	for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
 		auto ba_face = amrex::convert(ba, amrex::IntVect::TheDimensionVector(idim));
-		leftState[idim] = amrex::MultiFab(ba_face, dm, nvars, reconstructRange);
-		rightState[idim] = amrex::MultiFab(ba_face, dm, nvars, reconstructRange);
+		leftState[idim] = amrex::MultiFab(ba_face, dm, nvars, reconstructGhost);
+		rightState[idim] = amrex::MultiFab(ba_face, dm, nvars, reconstructGhost);
 		flux[idim] = amrex::MultiFab(ba_face, dm, nvars, 0);
 		facevel[idim] = amrex::MultiFab(ba_face, dm, 1, 0);
 	}
@@ -1153,15 +1207,50 @@ auto RadhydroSimulation<problem_t>::computeHydroFluxes(amrex::MultiFab const &co
 
 	// compute flux functions
 	AMREX_D_TERM(hydroFluxFunction<FluxDir::X1>(primVar, leftState[0], rightState[0], flux[0], facevel[0], flatCoefs[0], flatCoefs[1], flatCoefs[2],
-						    reconstructRange, nvars);
+						    reconstructGhost, nvars);
 		     , hydroFluxFunction<FluxDir::X2>(primVar, leftState[1], rightState[1], flux[1], facevel[1], flatCoefs[0], flatCoefs[1], flatCoefs[2],
-						      reconstructRange, nvars);
+						      reconstructGhost, nvars);
 		     , hydroFluxFunction<FluxDir::X3>(primVar, leftState[2], rightState[2], flux[2], facevel[2], flatCoefs[0], flatCoefs[1], flatCoefs[2],
-						      reconstructRange, nvars);)
+						      reconstructGhost, nvars);)
 
 	// synchronization point to prevent MultiFabs from going out of scope
 	amrex::Gpu::streamSynchronizeAll();
 
+	// LOW LEVEL DEBUGGING: output all of the temporary MultiFabs
+	if (lowLevelDebuggingOutput_ == 1) {
+		// write primitive cell-centered state
+		std::string plotfile_name = CustomPlotFileName("debug_reconstruction", istep[lev] + 1);
+		WriteSingleLevelPlotfile(plotfile_name, primVar, componentNames_cc_, geom[lev], 0.0, istep[lev] + 1);
+
+		// write flattening coefficients
+		std::string flatx_filename = CustomPlotFileName("debug_flattening_x", istep[lev] + 1);
+		std::string flaty_filename = CustomPlotFileName("debug_flattening_y", istep[lev] + 1);
+		std::string flatz_filename = CustomPlotFileName("debug_flattening_z", istep[lev] + 1);
+		amrex::Vector<std::string> flatCompNames{"chi"};
+		WriteSingleLevelPlotfile(flatx_filename, flatCoefs[0], flatCompNames, geom[lev], 0.0, istep[lev] + 1);
+		WriteSingleLevelPlotfile(flaty_filename, flatCoefs[1], flatCompNames, geom[lev], 0.0, istep[lev] + 1);
+		WriteSingleLevelPlotfile(flatz_filename, flatCoefs[2], flatCompNames, geom[lev], 0.0, istep[lev] + 1);
+
+		// write L interface states
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			if (amrex::ParallelDescriptor::IOProcessor()) {
+				std::filesystem::create_directories(plotfile_name + "/raw_fields/Level_" + std::to_string(lev));
+			}
+			std::string fullprefix =
+			    amrex::MultiFabFileFullPrefix(lev, plotfile_name, "raw_fields/Level_", std::string("StateL_") + quokka::face_dir_str[idim]);
+			amrex::VisMF::Write(leftState[idim], fullprefix);
+		}
+		// write R interface states
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			if (amrex::ParallelDescriptor::IOProcessor()) {
+				std::filesystem::create_directories(plotfile_name + "/raw_fields/Level_" + std::to_string(lev));
+			}
+			std::string fullprefix =
+			    amrex::MultiFabFileFullPrefix(lev, plotfile_name, "raw_fields/Level_", std::string("StateR_") + quokka::face_dir_str[idim]);
+			amrex::VisMF::Write(rightState[idim], fullprefix);
+		}
+	}
+
 	// return flux and face-centered velocities
 	return std::make_pair(std::move(flux), std::move(facevel));
 }

src/simulation.hpp

@@ -216,6 +216,8 @@ template <typename problem_t> class AMRSimulation : public amrex::AmrCore
 
 	// I/O functions
 	[[nodiscard]] auto PlotFileName(int lev) const -> std::string;
+	[[nodiscard]] auto CustomPlotFileName(const char *base, int lev) const -> std::string;
+	[[nodiscard]] auto GetPlotfileVarNames() const -> amrex::Vector<std::string>;
 	[[nodiscard]] auto PlotFileMF() const -> amrex::Vector<amrex::MultiFab>;
 	[[nodiscard]] auto PlotFileMFAtLevel(int lev) const -> amrex::MultiFab;
 	void WriteMetadataFile(std::string const &plotfilename) const;
@@ -1540,6 +1542,13 @@ template <typename problem_t> void AMRSimulation<problem_t>::AverageDownTo(int c
 // get plotfile name
 template <typename problem_t> auto AMRSimulation<problem_t>::PlotFileName(int lev) const -> std::string { return amrex::Concatenate(plot_file, lev, 5); }
 
+// get plotfile name
+template <typename problem_t> auto AMRSimulation<problem_t>::CustomPlotFileName(const char *base, int lev) const -> std::string
+{
+	std::string base_str(base);
+	return amrex::Concatenate(base_str, lev, 5);
+}
+
 template <typename problem_t>
 void AMRSimulation<problem_t>::AverageFCToCC(amrex::MultiFab &mf_cc, const amrex::MultiFab &mf_fc, int idim, int dstcomp_start, int srccomp_start,
 					     int srccomp_total, int nGrow) const
@@ -1690,6 +1699,19 @@ template <typename problem_t> void AMRSimulation<problem_t>::RenderAscent()
 }
 #endif // AMREX_USE_ASCENT
 
+template <typename problem_t> auto AMRSimulation<problem_t>::GetPlotfileVarNames() const -> amrex::Vector<std::string>
+{
+	amrex::Vector<std::string> varnames;
+	varnames.insert(varnames.end(), componentNames_cc_.begin(), componentNames_cc_.end());
+	if constexpr (Physics_Indices<problem_t>::nvarTotal_fc > 0) {
+		for (int icomp = 0; icomp < Physics_Indices<problem_t>::nvarTotal_fc; ++icomp) {
+			varnames.push_back(componentNames_fc_[icomp]);
+		}
+	}
+	varnames.insert(varnames.end(), derivedNames_.begin(), derivedNames_.end());
+	return varnames;
+}
+
 // write plotfile to disk
 template <typename problem_t> void AMRSimulation<problem_t>::WritePlotFile() const
 {
@@ -1703,17 +1725,11 @@ template <typename problem_t> void AMRSimulation<problem_t>::WritePlotFile() con
 #endif
 
 	// now construct output and submit to async write queue
-	const std::string &plotfilename = PlotFileName(istep[0]);
 	amrex::Vector<amrex::MultiFab> mf = PlotFileMF();
 	amrex::Vector<const amrex::MultiFab *> mf_ptr = amrex::GetVecOfConstPtrs(mf);
-	amrex::Vector<std::string> varnames;
-	varnames.insert(varnames.end(), componentNames_cc_.begin(), componentNames_cc_.end());
-	if constexpr (Physics_Indices<problem_t>::nvarTotal_fc > 0) {
-		for (int icomp = 0; icomp < Physics_Indices<problem_t>::nvarTotal_fc; ++icomp) {
-			varnames.push_back(componentNames_fc_[icomp]);
-		}
-	}
-	varnames.insert(varnames.end(), derivedNames_.begin(), derivedNames_.end());
+
+	const std::string &plotfilename = PlotFileName(istep[0]);
+	auto varnames = GetPlotfileVarNames();
 
 	// write plotfile
 	amrex::Print() << "Writing plotfile " << plotfilename << "\n";

tests/ShockCloud_32.in

@@ -0,0 +1,53 @@
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo     =  0.0  0.0  0.0 
+geometry.prob_hi     =  2.4688e+21  6.172e+20  6.172e+20  # 800 x 200 x 200 pc
+geometry.is_periodic =  0    1    1
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v               = 1     # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell          = 64 16 16
+amr.max_level       = 0     # number of levels = max_level + 1
+amr.blocking_factor_x = 64
+amr.blocking_factor_y = 16
+amr.blocking_factor_z = 16  # grid size must be divisible by this
+amr.max_grid_size     = 64
+
+# *****************************************************************
+# Quokka options
+# *****************************************************************
+cfl = 0.3
+do_reflux = 1
+do_subcycle = 1
+max_walltime = 6:00:00
+max_timesteps = 350
+checkpoint_interval = 5000
+plotfile_interval = 100
+ascent_interval = 1
+derived_vars = pressure entropy nH temperature cooling_length \
+	       cloud_fraction lab_velocity_x mass velocity_mag
+
+hydro.rk_integrator_order = 1
+hydro.reconstruction_order = 1
+hydro.use_dual_energy = 0
+hydro.low_level_debugging_output = 1
+
+cooling.enabled = 0
+cooling.read_tables_even_if_disabled = 1
+cooling.grackle_data_file = "../extern/grackle_data_files/input/CloudyData_UVB=HM2012.h5"
+temperature_floor = 100
+
+sharp_cloud_edge = 1
+do_frame_shift = 0
+nH_bg    = 0.003356403333 # cm^-3
+nH_cloud = 0.003356403333 # cm^-3
+P_over_k = 1.304005e+04   # K cm^-3
+R_cloud_pc = 16.09084149928867 # pc
+Mach_shock = 1.2