add TabulatedCooling test problem

src/ShockCloud/CMakeLists.txt

@@ -1,10 +1,12 @@
 if (AMReX_SPACEDIM GREATER_EQUAL 3)
     add_executable(shock_cloud cloud.cpp ${QuokkaObjSources})
     add_executable(test_grackle_cooling test_grackle_cooling.cpp ${QuokkaObjSources})
+    add_executable(test_cloudy_cooling test_cloudy_cooling.cpp ${QuokkaObjSources})
 
     if(AMReX_GPU_BACKEND MATCHES "CUDA")
         setup_target_for_cuda_compilation(shock_cloud)
         setup_target_for_cuda_compilation(test_grackle_cooling)
+        setup_target_for_cuda_compilation(test_cloudy_cooling)
     endif(AMReX_GPU_BACKEND MATCHES "CUDA")
 
 endif(AMReX_SPACEDIM GREATER_EQUAL 3)

src/ShockCloud/test_cloudy_cooling.cpp

@@ -0,0 +1,106 @@
+//==============================================================================
+// 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 cloud.cpp
+/// \brief Implements a shock-cloud problem with radiative cooling.
+///
+
+// uncomment this to debug the root-finding code (does NOT work on GPU!)
+// #define BOOST_MATH_INSTRUMENT
+
+#include "AMReX_ParmParse.H"
+#include "EOS.hpp"
+#include "ODEIntegrate.hpp"
+#include "TabulatedCooling.hpp"
+
+using amrex::Real;
+using namespace quokka::TabulatedCooling;
+
+struct ShockCloud {
+}; // dummy type to allow compile-type polymorphism via template specialization
+
+template <> struct quokka::EOS_Traits<ShockCloud> {
+	static constexpr double gamma = 5. / 3.; // default value
+};
+
+auto problem_main() -> int
+{
+	// Problem parameters
+	const Real rho = 2.27766918428822386e-22;  // g cm^-3;
+	const Real Eint = 1.11777608454088878e-11; // erg cm^-3
+	const Real dt = 1.92399749834457487e8;	   // s
+
+	// Read Cloudy tables
+	cloudy_tables cloudyTables;
+	std::string filename;
+	amrex::ParmParse const pp("cooling");
+	pp.query("hdf5_data_file", filename);
+	readCloudyData(filename, cloudyTables);
+	auto tables = cloudyTables.const_tables();
+
+	const Real T = ComputeTgasFromEgas(rho, Eint, quokka::EOS_Traits<ShockCloud>::gamma, tables);
+
+	const Real rhoH = rho * cloudy_H_mass_fraction;
+	const Real nH = rhoH / (C::m_p + C::m_e);
+	const Real log_nH = std::log10(nH);
+
+	const Real C = (quokka::EOS_Traits<ShockCloud>::gamma - 1.) * Eint / (C::k_B * (rho / (C::m_p + C::m_e)));
+	const Real mu = interpolate2d(log_nH, std::log10(T), tables.log_nH, tables.log_Tgas, tables.meanMolWeight);
+	const Real relerr = std::abs((C * mu - T) / T);
+
+	printf("\nrho = %.17e, Eint = %.17e, mu = %f, Tgas = %e, relerr = %e\n", rho, Eint, mu, T, relerr); // NOLINT
+
+	const Real reltol_floor = 0.01;
+	const Real rtol = 1.0e-4;	// not recommended to change this
+	constexpr Real T_floor = 100.0; // K
+	const Real gamma = quokka::EOS_Traits<ShockCloud>::gamma;
+
+	ODEUserData user_data{rho, gamma, tables};
+	quokka::valarray<Real, 1> y = {Eint};
+	quokka::valarray<Real, 1> const abstol = {reltol_floor * ComputeEgasFromTgas(rho, T_floor, gamma, tables)};
+
+	// do integration with RK2 (Heun's method)
+	int nsteps = 0;
+	rk_adaptive_integrate(user_rhs, 0, y, dt, &user_data, rtol, abstol, nsteps);
+
+	// check if integration failed
+	if (nsteps >= maxStepsODEIntegrate) {
+		Real const T = ComputeTgasFromEgas(rho, Eint, quokka::EOS_Traits<ShockCloud>::gamma, tables);
+		Real const Edot = cloudy_cooling_function(rho, T, tables);
+		Real const t_cool = Eint / Edot;
+		printf("max substeps exceeded! rho = %g, Eint = %g, T = %g, cooling " // NOLINT
+		       "time = %g\n",
+		       rho, Eint, T, t_cool);
+	} else {
+		Real const T = ComputeTgasFromEgas(rho, Eint, quokka::EOS_Traits<ShockCloud>::gamma, tables);
+		Real const Edot = cloudy_cooling_function(rho, T, tables);
+		Real const t_cool = Eint / Edot;
+		printf("success! rho = %g, Eint = %g, T = %g, cooling time = %g\n", rho, Eint, T, t_cool); // NOLINT
+	}
+
+#if 0
+  // compute Field length
+  auto lambda_F = [=] (Real nH0, Real T0) {
+    const Real rho0 = nH0 * hydrogen_mass_cgs_ / cloudy_H_mass_fraction; // g cm^-3
+    const Real Edot = cloudy_cooling_function(rho0, T0, tables);
+    const Real ln_L = 29.7 + std::log(std::pow(nH0, -1./2.) * (T0 / 1.0e6));
+    const Real conductivity = 1.84e-5 * std::pow(T0, 5./2.) / ln_L;
+    const Real l = std::sqrt( conductivity * T0 / std::abs(Edot) );
+    return std::make_pair(Edot, l);
+  };
+
+  auto [Edot0, l0] = lambda_F(1.0, 4.0e5);
+  amrex::Print() << "Edot(nH = 1.0, T = 4e5) = " << Edot0 << "\n";
+  amrex::Print() << "lambda_F = " << (l0 / 3.086e18) << " pc\n\n";
+
+  auto [Edot1, l1] = lambda_F(0.1, 4.0e5);
+  amrex::Print() << "Edot(nH = 0.1, T = 4e5) = " << Edot1 << "\n";
+  amrex::Print() << "lambda_F = " << (l1 / 3.086e18) << " pc\n\n";
+#endif
+
+	// Cleanup and exit
+	const int status = 0;
+	return status;
+}

tests/TabulatedCooling.in

@@ -0,0 +1,26 @@
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo     =  0.0  0.0  0.0 
+geometry.prob_hi     =  1.11096e18  4.44384e18  1.11096e18
+geometry.is_periodic =  1    0    1
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v               = 0     # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell          = 128 512 128
+amr.max_level       = 0     # number of levels = max_level + 1
+amr.blocking_factor = 128   # grid size must be divisible by this
+amr.max_grid_size   = 128
+
+do_reflux = 0
+do_subcycle = 0
+
+cooling.enabled = 1
+cooling.hdf5_data_file = "./isrf_1000Go_grains.h5"
+temperature_floor = 10  # K