output 2D slice plotfiles (#562)

### Description
Output 2D slices of the simulation as AMReX plotfiles. The
implementation is adaped from PeleLMeX
([documentation](https://amrex-combustion.github.io/PeleLMeX/manual/html/LMeXControls.html#run-time-diagnostics)).

It can be turned on for any variables (including derived variables) with
runtime parameters, e.g.:
```
quokka.diagnostics = slice_z           # Space-separated name(s) of diagnostics (arbitrary)
quokka.slice_z.type = DiagFramePlane   # Diagnostic type (others may be added in the future)
quokka.slice_z.file = slicez_plt       # Output file prefix (should end in "plt")
quokka.slice_z.normal = 2              # Plane normal (0 == x, 1 == y, 2 == z)
quokka.slice_z.center = 2.4688e20      # Coordinate in the normal direction
quokka.slice_z.int    = 10             # Output interval (in number of coarse steps)
quokka.slice_z.interpolation = Linear  # Interpolation type: Linear or Quadratic (default: Linear)

# The problem must output these derived variable(s)
derived_vars = temperature
# List of variables to include in output
quokka.slice_z.field_names = gasDensity gasInternalEnergy temperature
```

### Related issues
Closes #237.

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [ ] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/CMakeLists.txt

@@ -121,16 +121,15 @@ link_libraries(yaml-cpp::yaml-cpp)
 
 include(CTest)
 
-#create an object library for files that should be compiled with all test problems
-set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" 
-                      "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp" "${openPMDSources}"
-                      "${gamma_law_sources}")
-#we don't use it anymore because it gives issues on Cray systems
-#add_library(QuokkaObj OBJECT ${QuokkaObjSources} ${gamma_law_sources})
-#if(AMReX_GPU_BACKEND MATCHES "CUDA")
-#    setup_target_for_cuda_compilation(QuokkaObj)
-#endif()
-#link_libraries(QuokkaObj)
+set (QuokkaSourcesNoEOS "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp"
+                        "${CMAKE_CURRENT_SOURCE_DIR}/DiagBase.cpp" 
+                        "${CMAKE_CURRENT_SOURCE_DIR}/DiagFilter.cpp" 
+                        "${CMAKE_CURRENT_SOURCE_DIR}/DiagFramePlane.cpp" 
+                        "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" 
+                        "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp"
+                        "${openPMDSources}")
+
+set (QuokkaObjSources "${QuokkaSourcesNoEOS}" "${gamma_law_sources}")
 
 add_subdirectory(Advection)
 add_subdirectory(Advection2D)

src/DiagBase.H

@@ -0,0 +1,44 @@
+#ifndef DIAGBASE_H
+#define DIAGBASE_H
+
+#include "AMReX_MultiFab.H"
+#include "DiagFilter.H"
+#include "Factory.H"
+
+class DiagBase : public quokka::Factory<DiagBase>
+{
+      public:
+	~DiagBase() override = default;
+
+	static auto base_identifier() -> std::string { return "DiagBase"; }
+
+	virtual void init(const std::string &a_prefix, std::string_view a_diagName);
+
+	virtual void close() = 0;
+
+	auto needUpdate() const -> bool { return need_update; }
+
+	virtual auto doDiag(const amrex::Real &a_time, int a_nstep) -> bool;
+
+	virtual void prepare(int a_nlevels, const amrex::Vector<amrex::Geometry> &a_geoms, const amrex::Vector<amrex::BoxArray> &a_grids,
+			     const amrex::Vector<amrex::DistributionMapping> &a_dmap, const amrex::Vector<std::string> &a_varNames);
+
+	virtual void processDiag(int a_nstep, const amrex::Real &a_time, const amrex::Vector<const amrex::MultiFab *> &a_state,
+				 const amrex::Vector<std::string> &a_varNames) = 0;
+
+	virtual void addVars(amrex::Vector<std::string> &a_varList);
+
+	static auto getFieldIndex(const std::string &a_field, const amrex::Vector<std::string> &a_varList) -> int;
+
+      protected:
+	std::string m_diagfile;
+	int m_verbose{0};
+	amrex::Real m_per{-1.0};
+	int m_interval{-1};
+	bool need_update{true};
+	bool first_time{true};
+	amrex::Vector<DiagFilter> m_filters{};
+	amrex::Gpu::DeviceVector<DiagFilterData> m_filterData;
+};
+
+#endif

src/DiagBase.cpp

@@ -0,0 +1,81 @@
+#include "DiagBase.H"
+#include "AMReX_ParmParse.H"
+
+void DiagBase::init(const std::string &a_prefix, std::string_view a_diagName)
+{
+	amrex::ParmParse const pp(a_prefix);
+
+	// IO
+	pp.query("int", m_interval);
+	pp.query("per", m_per);
+	m_diagfile = a_diagName;
+	pp.query("file", m_diagfile);
+	AMREX_ASSERT(m_interval > 0 || m_per > 0.0);
+
+	// Filters
+	int const nFilters = pp.countval("filters");
+	amrex::Vector<std::string> filtersName;
+	if (nFilters > 0) {
+		m_filters.resize(nFilters);
+		filtersName.resize(nFilters);
+	}
+	for (int n = 0; n < nFilters; ++n) {
+		pp.get("filters", filtersName[n], n);
+		const std::string filter_prefix = a_prefix + "." + filtersName[n];
+		m_filters[n].init(filter_prefix);
+	}
+}
+
+void DiagBase::prepare(int /*a_nlevels*/, const amrex::Vector<amrex::Geometry> & /*a_geoms*/, const amrex::Vector<amrex::BoxArray> & /*a_grids*/,
+		       const amrex::Vector<amrex::DistributionMapping> & /*a_dmap*/, const amrex::Vector<std::string> &a_varNames)
+{
+	if (first_time) {
+		int const nFilters = static_cast<int>(m_filters.size());
+		// Move the filter data to the device
+		for (int n = 0; n < nFilters; ++n) {
+			m_filters[n].setup(a_varNames);
+		}
+		amrex::Vector<DiagFilterData> hostFilterData;
+		for (int n = 0; n < nFilters; ++n) {
+			hostFilterData.push_back(m_filters[n].m_fdata);
+		}
+		m_filterData.resize(nFilters);
+		amrex::Gpu::copy(amrex::Gpu::hostToDevice, hostFilterData.begin(), hostFilterData.end(), m_filterData.begin());
+	}
+}
+
+auto DiagBase::doDiag(const amrex::Real &a_time, int a_nstep) -> bool
+{
+	bool willDo = false;
+	if (m_interval > 0 && (a_nstep % m_interval == 0)) {
+		willDo = true;
+	}
+
+	// TODO(bwibking): output based on a_time
+	amrex::ignore_unused(a_time);
+
+	return willDo;
+}
+
+void DiagBase::addVars(amrex::Vector<std::string> &a_varList)
+{
+	int const nFilters = static_cast<int>(m_filters.size());
+	for (int n = 0; n < nFilters; ++n) {
+		a_varList.push_back(m_filters[n].m_filterVar);
+	}
+}
+
+auto DiagBase::getFieldIndex(const std::string &a_field, const amrex::Vector<std::string> &a_varList) -> int
+{
+	int index = -1;
+	for (int n{0}; n < a_varList.size(); ++n) {
+		if (a_varList[n] == a_field) {
+			index = n;
+			break;
+		}
+	}
+	if (index < 0) {
+		amrex::Abort("Field : " + a_field + " wasn't found in available fields");
+	}
+	return index;
+}

src/DiagFilter.H

@@ -0,0 +1,21 @@
+#ifndef DIAGFILTERS_H
+#define DIAGFILTERS_H
+
+#include <AMReX.H>
+#include <AMReX_Utility.H>
+
+struct DiagFilterData {
+	int m_filterVarIdx{-1};
+	amrex::Real m_low_val{0.0};
+	amrex::Real m_high_val{0.0};
+};
+
+struct DiagFilter {
+	std::string m_filterVar;
+
+	void init(const std::string &a_prefix);
+	void setup(const amrex::Vector<std::string> &a_varNames);
+
+	DiagFilterData m_fdata;
+};
+#endif

src/DiagFilter.cpp

@@ -0,0 +1,45 @@
+#include "DiagFilter.H"
+#include "AMReX_ParmParse.H"
+
+void DiagFilter::init(const std::string &a_prefix)
+{
+	amrex::ParmParse const pp(a_prefix);
+
+	pp.query("field_name", m_filterVar);
+	if (m_filterVar.empty()) {
+		amrex::Abort("filter: " + a_prefix + " is missing a field_name !");
+	}
+
+	int definedRange = 0;
+	if (pp.countval("value_greater") != 0) {
+		pp.get("value_greater", m_fdata.m_low_val);
+		m_fdata.m_high_val = AMREX_REAL_MAX;
+		definedRange = 1;
+	} else if (pp.countval("value_less") != 0) {
+		pp.get("value_less", m_fdata.m_high_val);
+		m_fdata.m_low_val = AMREX_REAL_LOWEST;
+		definedRange = 1;
+	} else if (pp.countval("value_inrange") != 0) {
+		amrex::Vector<amrex::Real> range{0.0};
+		pp.getarr("value_inrange", range, 0, 2);
+		m_fdata.m_low_val = std::min(range[0], range[1]);
+		m_fdata.m_high_val = std::max(range[0], range[1]);
+		definedRange = 1;
+	}
+	if (definedRange == 0) {
+		amrex::Abort("filter: " + a_prefix + " is missing a range definition !");
+	}
+}
+
+void DiagFilter::setup(const amrex::Vector<std::string> &a_varNames)
+{
+	m_fdata.m_filterVarIdx = -1;
+	for (int n{0}; n < a_varNames.size(); ++n) {
+		if (a_varNames[n] == m_filterVar) {
+			m_fdata.m_filterVarIdx = n;
+		}
+	}
+	if (m_fdata.m_filterVarIdx < 0) {
+		amrex::Abort("Couldn't find filter field: " + m_filterVar);
+	}
+}

src/DiagFramePlane.H

@@ -0,0 +1,66 @@
+#ifndef DIAGFRAMEPLANE_H
+#define DIAGFRAMEPLANE_H
+
+#include "AMReX_VisMF.H"
+#include "DiagBase.H"
+
+class DiagFramePlane : public DiagBase::Register<DiagFramePlane>
+{
+      public:
+	static auto identifier() -> std::string { return "DiagFramePlane"; }
+
+	enum InterpType { Linear, Quadratic };
+
+	void init(const std::string &a_prefix, std::string_view a_diagName) override;
+
+	void prepare(int a_nlevels, const amrex::Vector<amrex::Geometry> &a_geoms, const amrex::Vector<amrex::BoxArray> &a_grids,
+		     const amrex::Vector<amrex::DistributionMapping> &a_dmap, const amrex::Vector<std::string> &a_varNames) override;
+
+	void processDiag(int a_nstep, const amrex::Real &a_time, const amrex::Vector<const amrex::MultiFab *> &a_state,
+			 const amrex::Vector<std::string> &a_varNames) override;
+
+	void addVars(amrex::Vector<std::string> &a_varList) override;
+
+	void Write2DMultiLevelPlotfile(const std::string &a_pltfile, int a_nlevels, const amrex::Vector<const amrex::MultiFab *> &a_slice,
+				       const amrex::Vector<std::string> &a_varnames, const amrex::Vector<amrex::Geometry> &a_geoms, const amrex::Real &a_time,
+				       const amrex::Vector<int> &a_steps, const amrex::Vector<amrex::IntVect> &a_rref);
+
+	static void VisMF2D(const amrex::MultiFab &a_mf, const std::string &a_mf_name);
+
+	static void Write2DMFHeader(const std::string &a_mf_name, amrex::VisMF::Header &hdr, int coordinatorProc, MPI_Comm comm);
+
+	static void Find2FOffsets(const amrex::FabArray<amrex::FArrayBox> &mf, const std::string &filePrefix, amrex::VisMF::Header &hdr,
+				  amrex::VisMF::Header::Version /*whichVersion*/, amrex::NFilesIter &nfi, int nOutFile, MPI_Comm comm);
+
+	static void write_2D_header(std::ostream &os, const amrex::FArrayBox &f, int nvar);
+
+	static void Write2DPlotfileHeader(std::ostream &HeaderFile, int nlevels, const amrex::Vector<amrex::BoxArray> &bArray,
+					  const amrex::Vector<std::string> &varnames, const amrex::Vector<amrex::Geometry> &geom, const amrex::Real &time,
+					  const amrex::Vector<int> &level_steps, const amrex::Vector<amrex::IntVect> &ref_ratio,
+					  const std::string &versionName = "HyperCLaw-V1.1", const std::string &levelPrefix = "Level_",
+					  const std::string &mfPrefix = "Cell");
+
+	void close() override {}
+
+      private:
+	// Variables output
+	amrex::Vector<std::string> m_fieldNames;
+	amrex::Gpu::DeviceVector<int> m_fieldIndices_d;
+
+	// Plane definition
+	int m_normal;
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> m_center;
+
+	// Interpolation data
+	InterpType m_interpType;
+	amrex::Vector<amrex::GpuArray<amrex::Real, 3>> m_intwgt;
+	amrex::Vector<int> m_k0;
+
+	// 2D-plane boxArray vector
+	amrex::Geometry m_geomLev0;
+	amrex::Vector<amrex::BoxArray> m_sliceBA;
+	amrex::Vector<amrex::Vector<int>> m_dmConvert;
+	amrex::Vector<amrex::DistributionMapping> m_sliceDM;
+};
+
+#endif