diff --git a/src/ShockCloud/CMakeLists.txt b/src/ShockCloud/CMakeLists.txt index 5d3bb12d3..4f4cc286f 100644 --- a/src/ShockCloud/CMakeLists.txt +++ b/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) diff --git a/src/ShockCloud/test_cloudy_cooling.cpp b/src/ShockCloud/test_cloudy_cooling.cpp new file mode 100644 index 000000000..c624fa5ec --- /dev/null +++ b/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 { + 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::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::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::gamma; + + ODEUserData user_data{rho, gamma, tables}; + quokka::valarray y = {Eint}; + quokka::valarray 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::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::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; +} diff --git a/tests/TabulatedCooling.in b/tests/TabulatedCooling.in new file mode 100644 index 000000000..c6d2189a1 --- /dev/null +++ b/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