diff --git a/.gitmodules b/.gitmodules index 4788484d..c1eebde1 100644 --- a/.gitmodules +++ b/.gitmodules @@ -16,3 +16,6 @@ [submodule "external/Catch2"] path = external/Catch2 url = https://github.com/catchorg/Catch2.git +[submodule "external/kokkos-kernels"] + path = external/kokkos-kernels + url = https://github.com/kokkos/kokkos-kernels diff --git a/CMakeLists.txt b/CMakeLists.txt index 8efaf728..2d2f0427 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -176,6 +176,10 @@ set(CMAKE_FIND_PACKAGE_NO_SYSTEM_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE) message("Configuring singularity-eos") add_subdirectory(external/singularity-eos singularity-eos) +# Kokkos kernels +message("Configuring Kokkos Kernels") +add_subdirectory(external/kokkos-kernels) + # REBOUND message("Configuring REBOUND") # Define the source and destination directories diff --git a/external/kokkos-kernels b/external/kokkos-kernels new file mode 160000 index 00000000..2899ebab --- /dev/null +++ b/external/kokkos-kernels @@ -0,0 +1 @@ +Subproject commit 2899ebaba1849b47819ffec5d30672119ace7712 diff --git a/inputs/atmosphere/escape_1d.in b/inputs/atmosphere/escape_1d.in new file mode 100644 index 00000000..99a935d2 --- /dev/null +++ b/inputs/atmosphere/escape_1d.in @@ -0,0 +1,95 @@ +# ============================================================================= +# Artemis input deck: 1D spherical planetary atmospheric escape +# +# Two species (H, mu=1 and O, mu=16) in a 1D spherically-symmetric atmosphere +# under point-mass gravity. Chapman-Cowling (hard-sphere) drag couples them. +# Species 0 (H) is light enough to escape; species 1 (O) is gravitationally +# bound but gets entrained by species 0 through drag. + + + problem = escape_1d + coordinates = spherical + radial_spacing = logarithmic + + + problem_id = drag + + + nlim = 1 + tlim = 50.0 # ~14 t_dyn (t_dyn = r0/v_esc = 1/sqrt(2)) + integrator = rk2 + ncycle_out = 1 + + + nx1 = 64 + x1min = 0.0 # inner radius r0 + x1max = 2.302585092994046 # outer radius + ix1_bc = hydrostatic # custom inner BC (hydrostatic inflow) + ox1_bc = outflow # custom outer BC (outflow) + + nx2 = 1 + x2min = 0.0 + x2max = 3.14159265358979 + ix2_bc = outflow + ox2_bc = outflow + + nx3 = 1 + x3min = -0.5 + x3max = 0.5 + ix3_bc = outflow + ox3_bc = outflow + + + nx1 = 64 + nx2 = 1 + nx3 = 1 + + + file_type = hdf5 + dt = 0.1 + variables = gas.prim.density, gas.prim.velocity, gas.prim.sie, gas.prim.temperature, gas.prim.pressure + + + file_type = hst + dt = 0.01 + +# Gas: two-species ideal gas + + nspecies = 2 + gamma = 1.6666666667, 1.6666666667 + mu = 1.0, 16.0 # H, O [AMU] + dfloor = 1.0e-12 + siefloor = 1.0e-12 + +# Gravity: point mass at origin + + mass = 1.0 # G*M in code units + + + + gas = true + gravity = true + drag = true + +# Drag: full Chapman-Cowling coupling + + type = full + + + collision_model = hard_sphere + # mu [AMU], sigma [code length = r0], dof (H/O treated here as monatomic) + mu = 1.0 , 16.0 + sigma = 1.0e-3, 1.0e-3 # collision diameters ~1e-3 r0 (calibrate for your scenario) + dof = 3.0, 3.0 + +# Problem parameters + + # Isothermal temperature in EOS/code temperature units for this scalefree deck. + # T0 = 0.2 yields lambda_H = mu_H/T0 = 5. The inner BC maintains hydrostatic + # density/SIE and extrapolates outward radial velocity; it does not impose a + # Parker-wind velocity at r0. + T0 = 0.2 # code units (NOT Kelvin) + r0 = 1.0 # inner radius [code length] + gm = 1.0 # G*M [code units] + rho0_0 = 1.0e-4 # base density species 0 (H) at r0 + rho0_1 = 1.0e-3 # base density species 1 (O) at r0 (heavier, higher abundance) diff --git a/inputs/atmosphere/escape_n3_1d.in b/inputs/atmosphere/escape_n3_1d.in new file mode 100644 index 00000000..7cde1291 --- /dev/null +++ b/inputs/atmosphere/escape_n3_1d.in @@ -0,0 +1,106 @@ +# ============================================================================= +# Artemis input deck: 1D spherical planetary atmospheric escape +# +# Two species (H, mu=1 and O, mu=16) in a 1D spherically-symmetric atmosphere +# under point-mass gravity. Chapman-Cowling (hard-sphere) drag couples them. +# Species 0 (H) is light enough to escape; species 1 (O) is gravitationally +# bound but gets entrained by species 0 through drag. + + + problem = escape_1d + coordinates = spherical + radial_spacing = logarithmic + + + enable_sparse = true + dealloc_count = 1 + alloc_threshold = 1e-8 + dealloc_threshold = 1e-10 + + + problem_id = drag + + + nlim = 100 + tlim = 50.0 # ~14 t_dyn (t_dyn = r0/v_esc = 1/sqrt(2)) + integrator = rk2 + ncycle_out = 1 + + + nx1 = 512 + x1min = 0.0 # inner radius r0 + x1max = 2.302585092994046 # outer radius + ix1_bc = hydrostatic # custom inner BC (hydrostatic inflow) + ox1_bc = outflow # custom outer BC (outflow) + + nx2 = 1 + x2min = 0.0 + x2max = 3.14159265358979 + ix2_bc = outflow + ox2_bc = outflow + + nx3 = 1 + x3min = -0.5 + x3max = 0.5 + ix3_bc = outflow + ox3_bc = outflow + + + nx1 = 8 + nx2 = 1 + nx3 = 1 + + + file_type = hdf5 +# dt = 0.1 + dn = 1 + variables = gas.prim.density, gas.prim.velocity, gas.prim.sie, gas.prim.temperature, gas.prim.pressure + use_final_label = no + + + file_type = hst + dt = 0.01 + + + nspecies = 3 + gamma = 1.6666666667, 1.6666666667, 1.6666666667 + mu = 1.0, 4.0026, 16.0 # H, He, O [AMU] + dfloor = 1.0e-12 + siefloor = 1.0e-12 + scr_level = 1 + density_allocation_threshold = 1e-8 + density_deallocation_threshold = 1e-10 + +# Gravity: point mass at origin + + mass = 1.0 # G*M in code units + + + + gas = true + gravity = true + drag = true + +# Drag: full Chapman-Cowling coupling + + type = full + + + collision_model = hard_sphere + # mu [AMU], sigma [code length = r0], dof (H/O treated here as monatomic) + mu = 1.0 , 4.0026, 16.0 + sigma = 1.0e-3, 1e-3, 1.0e-3 # collision diameters ~1e-3 r0 (calibrate for your scenario) + dof = 3.0, 3.0, 3.0 + +# Problem parameters + + # Isothermal temperature in EOS/code temperature units for this scalefree deck. + # T0 = 0.2 yields lambda_H = mu_H/T0 = 5. The inner BC maintains hydrostatic + # density/SIE and extrapolates outward radial velocity; it does not impose a + # Parker-wind velocity at r0. + T0 = 0.2 # code units (NOT Kelvin) + r0 = 1.0 # inner radius [code length] + gm = 1.0 # G*M [code units] + rho0_0 = 1.0e-4 # base density species 0 (H) at r0 + rho0_1 = 3.0e-4 # base density species 1 (He) at r0 + rho0_2 = 1.0e-3 # base density species 1 (O) at r0 (heavier, higher abundance) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 51386b94..df678137 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,6 +23,8 @@ set (SRC_LIST drag/drag.cpp drag/drag.hpp + drag/drag_impl.hpp + drag/collision_integrals.hpp dust/dust.cpp dust/dust.hpp @@ -137,10 +139,6 @@ if (ARTEMIS_PORTABLE_RESTART) add_compile_definitions(PORTABLE_RESTART) endif() -if (ARTEMIS_ENABLE_SESAME_EOS) - add_compile_definitions(WITH_SESAME) -endif() - # Generate library add_library(artemislib ${SRC_LIST}) @@ -166,6 +164,7 @@ endif() # Link libraries target_link_libraries(artemislib PUBLIC parthenon + KokkosKernels::kokkoskernels singularity-eos ${CMAKE_BINARY_DIR}/rebound/librebound.so lib_jaybenne_package diff --git a/src/artemis.cpp b/src/artemis.cpp index 7139c290..cf409b36 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -102,6 +102,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { const bool do_radiation = pin->GetOrAddBoolean("physics", "radiation", false); const bool do_coagulation = pin->GetOrAddBoolean("physics", "coagulation", false); const bool do_raytrace = pin->GetOrAddBoolean("physics", "raytrace", false); + const bool do_closure = pin->GetOrAddBoolean("physics", "closure", false); const bool do_orbital_advection = pin->GetOrAddBoolean("physics", "orbital_advection", false); @@ -145,6 +146,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { artemis->AddParam("do_orbital_advection", do_orbital_advection); artemis->AddParam("do_raytrace", do_raytrace); artemis->AddParam("update_fluxes", update_fluxes); + artemis->AddParam("do_closure", do_closure); // Set coordinate system const int ndim = ProblemDimension(pin.get()); @@ -172,7 +174,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_rotating_frame || do_orbital_advection) packages.Add(RotatingFrame::Initialize(pin.get())); if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get())); - if (do_drag) packages.Add(Drag::Initialize(pin.get())); + if (do_drag) packages.Add(Drag::Initialize(pin.get(), constants, packages)); // Operator split dust coagulation if (do_coagulation) { @@ -188,10 +190,11 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { packages.Add(Radiation::Initialize(pin.get(), constants, do_imc)); // Select between Jaybenne IMC or Moments if (do_imc) { - auto eos_h = packages.Get("gas")->Param("eos_h"); + auto eos_h_arr = + packages.Get("gas")->Param>("eos_h"); auto opacity_h = packages.Get("gas")->Param("opacity_h"); auto scattering_h = packages.Get("gas")->Param("scattering_h"); - packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, + packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h_arr(0), "radiation/imc")); PARTHENON_REQUIRE(coords == Coordinates::cartesian, "Jaybenne currently supports only Cartesian coordinates!"); @@ -238,6 +241,8 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_gas) Gas::AddHistory(coords, packages.Get("gas")->AllParams()); if (do_dust) Dust::AddHistory(coords, packages.Get("dust")->AllParams()); + params.Add("do_sparse", pin->GetBoolean("parthenon/sparse", "enable_sparse")); + // Add artemis package packages.Add(artemis); diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a4af03db..6c2f5d50 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -75,11 +75,13 @@ ArtemisDriver::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in do_moment = artemis_pkg->template Param("do_moment"); do_coagulation = artemis_pkg->template Param("do_coagulation"); do_raytrace = artemis_pkg->template Param("do_raytrace"); + do_closure = artemis_pkg->template Param("do_closure"); // Update fluxes option--gas fields are needed for radiation temperature updates but for // rad-only test problems turn off advection update_fluxes = artemis_pkg->template Param("update_fluxes"); + do_sparse = artemis_pkg->template Param("do_sparse"); ndim = pm->ndim; // Moments integrator @@ -384,6 +386,10 @@ TaskCollection ArtemisDriver::PostStepTasks() { if (pmesh->adaptive) { refine = tl.AddTask(new_dt, parthenon::Refinement::Tag>, u0.get()); } + auto dealloc = refine; + if (do_sparse) { + dealloc = tl.AddTask(refine, parthenon::Update::SparseDealloc, u0.get()); + } } return tc; diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index 0ac79ce4..72300c63 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -57,8 +57,9 @@ class ArtemisDriver : public EvolutionDriver { bool do_gravity, do_self_gravity, do_nbody, do_rotating_frame; bool do_orbital_advection; bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion; - bool do_coagulation, do_raytrace; + bool do_coagulation, do_raytrace, do_closure; bool update_fluxes; + bool do_sparse; const bool is_restart; }; diff --git a/src/derived/fill_derived.cpp b/src/derived/fill_derived.cpp index b65d13e3..092cb251 100644 --- a/src/derived/fill_derived.cpp +++ b/src/derived/fill_derived.cpp @@ -240,12 +240,12 @@ void PrimToCons(T *md) { // Extract gas parameters Real dflr_gas = Null(); Real sieflr_gas = Null(); - EOS eos_d; + ParArray1D eos_d; if (do_gas) { auto &gas_pkg = pm->packages.Get("gas"); dflr_gas = gas_pkg->template Param("dfloor"); sieflr_gas = gas_pkg->template Param("siefloor"); - eos_d = gas_pkg->template Param("eos_d"); + eos_d = gas_pkg->template Param>("eos_d"); } // Extract dust parameters @@ -291,7 +291,6 @@ void PrimToCons(T *md) { const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); if (do_gas) { - Real lambda[ArtemisUtils::lambda_max_vals] = {Null()}; for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { // Sync conserved and primitive density Real &w_d = vmesh(b, gas::prim::density(n), k, j, i); @@ -320,9 +319,9 @@ void PrimToCons(T *md) { const bool siefloor = (w_s > sieflr_gas); w_s = (siefloor)*w_s + (!siefloor) * sieflr_gas; u_u = w_s * u_d; - w_p = eos_d.PressureFromDensityInternalEnergy(w_d, w_s, lambda); - w_b = eos_d.BulkModulusFromDensityInternalEnergy(w_d, w_s, lambda); - w_t = eos_d.TemperatureFromDensityInternalEnergy(w_d, w_s, lambda); + w_p = eos_d(n).PressureFromDensityInternalEnergy(w_d, w_s); + w_b = eos_d(n).BulkModulusFromDensityInternalEnergy(w_d, w_s); + w_t = eos_d(n).TemperatureFromDensityInternalEnergy(w_d, w_s); // Sync conserved total energy const Real ke = 0.5 * w_d * (SQR(vel1) + SQR(vel2) + SQR(vel3)); diff --git a/src/drag/collision_integrals.hpp b/src/drag/collision_integrals.hpp new file mode 100644 index 00000000..ec9df68f --- /dev/null +++ b/src/drag/collision_integrals.hpp @@ -0,0 +1,136 @@ +//======================================================================================== +// (C) (or copyright) 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +#ifndef DRAG_COLLISION_INTEGRALS_HPP_ +#define DRAG_COLLISION_INTEGRALS_HPP_ + +//! \file collision_integrals.hpp +//! \brief Chapman-Cowling collision-integral backend for multifluid momentum and +//! energy coupling. These use empirical fits from Neufeld et al. 1972 and Kim & +//! Monroe 2014 + +#include + +#include "artemis.hpp" +#include "drag.hpp" + +namespace Drag { +namespace CollisionIntegrals { + +// --------------------------------------------------------------------------- +// Neufeld et al. (1972) polynomial fit coefficients for Omega^{(1,1)*} +// Valid for 0.3 <= T* <= 400. +// --------------------------------------------------------------------------- +namespace detail { +KOKKOS_FORCEINLINE_FUNCTION +Real Omega11Star(const Real Tstar) { + // Neufeld 1972 fit: A / T*^B + C/exp(D*T*) + E/exp(F*T*) + G/exp(H*T*) + constexpr Real A = 1.06036, B = 0.15610; + constexpr Real C = 0.19300, D = 0.47635; + constexpr Real E = 1.03587, F = 1.52996; + constexpr Real G = 1.76474, H = 3.89411; + return A / std::pow(Tstar, B) + C / std::exp(D * Tstar) + E / std::exp(F * Tstar) + + G / std::exp(H * Tstar); +} + +KOKKOS_FORCEINLINE_FUNCTION +Real Omega22Star(const Real Tstar) { + // Neufeld 1972 fit for Omega^{(2,2)*} + constexpr Real A = 1.16145, B = 0.14874; + constexpr Real C = 0.52487, D = 0.77320; + constexpr Real E = 2.16178, F = 2.43787; + constexpr Real G = -6.435e-4, H = 7.27371; + constexpr Real R = 18.0323, S = -0.76830, W = 7.27371; + // Extended fit (7-term, Kim & Monroe 2014 correction) + return A / std::pow(Tstar, B) + C / std::exp(D * Tstar) + E / std::exp(F * Tstar) + + G * std::sin(R * Tstar - W) / std::pow(Tstar, S); +} + +KOKKOS_FORCEINLINE_FUNCTION +Real sigma_ij(const Real sigma_i, const Real sigma_j) { + return 0.5 * (sigma_i + sigma_j); +} + +KOKKOS_FORCEINLINE_FUNCTION +Real eps_ij(const Real eps_i, const Real eps_j) { + return std::sqrt(std::max(eps_i * eps_j, 0.0)); +} + +KOKKOS_FORCEINLINE_FUNCTION +Real reduc_ij(const Real m_i, const Real m_j) { return m_i * m_j / (m_i + m_j); } + +} // namespace detail + +// --------------------------------------------------------------------------- +//! \fn DragCoeff +//! \brief Momentum-transfer rate coefficient K_{ij} [mass/(vol*time)] +//! \return K_{ij} such that d(rho_i v_i)/dt = -K_{ij}*(v_i - v_j) +// --------------------------------------------------------------------------- +KOKKOS_FORCEINLINE_FUNCTION +Real DragCoeff(const GasDragModel model, const Real kb_code, const Real m_i, + const Real m_j, const Real sig_i, const Real sig_j, const Real eps_i, + const Real eps_j, const Real rho_i, const Real rho_j, const Real T) { + using namespace detail; + constexpr Real pi = M_PI; + + const Real sij = sigma_ij(sig_i, sig_j); + Real O11 = 1.0; + if (model == GasDragModel::lj) { + const Real eij = eps_ij(eps_i, eps_j); + const Real Tstar = (eij > 0.0) ? T / eij : 1.0; + O11 = Omega11Star(std::max(Tstar, 0.3)); + } + + const Real v_rel = std::sqrt(8.0 * kb_code * T / (M_PI * reduc_ij(m_i, m_j))); + const Real K = + (4.0 / 3.0) * rho_i * rho_j / (m_i + m_j) * M_PI * SQR(sij) * v_rel * O11; + return std::max(K, 0.0); +} + +// --------------------------------------------------------------------------- +//! \fn ThermalRelaxRate +//! \brief Thermal energy exchange coefficient H_{ij} [energy/(vol*time*temperature)] +//! +//! Thermal exchange between species i and j: +//! dE_i/dt = H_{ij} * (T_j - T_i) +// --------------------------------------------------------------------------- +KOKKOS_FORCEINLINE_FUNCTION +Real ThermalRelaxRate(const GasDragModel model, const Real kb_code, const Real m_i, + const Real m_j, const Real sig_i, const Real sig_j, + const Real eps_i, const Real eps_j, const Real dof_i, + const Real dof_j, const Real cv_i, const Real cv_j, + const Real rho_i, const Real rho_j, const Real T) { + using namespace detail; + + const Real m_sum = m_i + m_j; + const Real mij = m_i * m_j / m_sum; + const Real K = + DragCoeff(model, kb_code, m_i, m_j, sig_i, sig_j, eps_i, eps_j, rho_i, rho_j, T); + + Real omega_E = 1.0; + if (model == GasDragModel::lj) { + const Real eij = eps_ij(eps_i, eps_j); + const Real Tstar = (eij > 0.0) ? T / eij : 1.0; + omega_E = Omega22Star(std::max(Tstar, 0.3)) / Omega11Star(std::max(Tstar, 0.3)); + } + + const Real chi = + (2.0 * mij / m_sum) * (2.0 * std::sqrt(dof_i * dof_j) / (dof_i + dof_j)); + const Real cv_pair = 2.0 * reduc_ij(cv_i, cv_j); + const Real H = K * omega_E * chi * cv_pair; + return std::max(H, 0.0); +} + +} // namespace CollisionIntegrals +} // namespace Drag + +#endif // DRAG_COLLISION_INTEGRALS_HPP_ diff --git a/src/drag/drag.cpp b/src/drag/drag.cpp index c1ce0255..c412cee0 100644 --- a/src/drag/drag.cpp +++ b/src/drag/drag.cpp @@ -14,15 +14,19 @@ // Artemis includes #include "drag.hpp" #include "artemis.hpp" +#include "drag_impl.hpp" #include "geometry/geometry.hpp" #include "utils/eos/eos.hpp" using ArtemisUtils::EOS; namespace Drag { + //---------------------------------------------------------------------------------------- //! \fn StateDescriptor Drag::Initialize //! \brief Adds intialization function for damping package -std::shared_ptr Initialize(ParameterInput *pin) { +std::shared_ptr Initialize(ParameterInput *pin, + const ArtemisUtils::Constants &constants, + Packages_t &packages) { auto drag = std::make_shared("drag"); Params ¶ms = drag->AllParams(); @@ -81,6 +85,17 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("stopping_time_params", StoppingTimeParams("dust/stopping_time", pin)); } + // Coupling::full — Chapman-Cowling N-fluid coupling + if (type == Coupling::full) { + PARTHENON_REQUIRE(do_gas, "drag type full requires at least one gas species"); + PARTHENON_REQUIRE(pin->GetOrAddInteger("gas", "nspecies", 1) > 1, + "drag type full requires gas.nspecies > 1"); + PARTHENON_REQUIRE(pin->DoesBlockExist("drag/full"), + "drag type full requires a [drag/full] input block"); + const auto gas_mu = packages.Get("gas")->Param>("mu"); + params.Add("full_coupling_params", FullCouplingParams(pin, constants, gas_mu)); + } + return drag; } @@ -96,7 +111,7 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { // Extract artemis parameters auto &artemis_pkg = pm->packages.Get("artemis"); const bool do_gas = artemis_pkg->template Param("do_gas"); - const bool do_dust = artemis_pkg->template Param("do_gas"); + const bool do_dust = artemis_pkg->template Param("do_dust"); // Extract drag parameters auto &drag_pkg = pm->packages.Get("drag"); @@ -109,13 +124,13 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { if (do_gas && gas_self_par.damp_to_visc) { auto &gas_pkg = pm->packages.Get("gas"); const auto &dp = gas_pkg->template Param("visc_params"); - const auto &eos_d = gas_pkg->template Param("eos_d"); + const auto eos_d = gas_pkg->template Param>("eos_d"); if (dp.type == Diffusion::DiffType::viscosity_plaw) { return SelfDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par); } else if (dp.type == Diffusion::DiffType::viscosity_alpha) { return SelfDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par); } else { PARTHENON_FAIL("The chosen viscosity model does not work with damping"); } @@ -128,7 +143,7 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { } else if (ctype == Coupling::simple_dust) { auto &gas_pkg = pm->packages.Get("gas"); auto &dust_pkg = pm->packages.Get("dust"); - const auto &eos_d = gas_pkg->template Param("eos_d"); + const auto eos_d = gas_pkg->template Param>("eos_d"); const auto stop_par = drag_pkg->template Param("stopping_time_params"); if (gas_self_par.damp_to_visc) { @@ -137,21 +152,21 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { if (stop_par.model == DragModel::constant) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } } else if (dp.type == Diffusion::DiffType::viscosity_alpha) { if (stop_par.model == DragModel::constant) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } } else { PARTHENON_FAIL("The chosen viscosity model does not work with damping"); @@ -160,18 +175,55 @@ TaskStatus DragSource(MeshData *md, const Real time, const Real dt) { Diffusion::DiffCoeffParams dp; if (stop_par.model == DragModel::constant) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } else if (stop_par.model == DragModel::stokes) { return SimpleDragSourceImpl( - md, time, dt, dp, eos_d, gas_self_par, dust_self_par, stop_par); + md, time, dt, dp, eos_d(0), gas_self_par, dust_self_par, stop_par); } } + } else if (ctype == Coupling::full) { + // Chapman-Cowling implicit coupling for N gas species. + // ApplyClosure assembles and solves the per-cell coupling matrix. + return ApplyClosure(md, dt); } else { PARTHENON_FAIL("Invalid drag model!"); + return TaskStatus::complete; } return TaskStatus::complete; } +template +TaskStatus ApplyClosure(MeshData *md, const Real dt) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Packing and indexing + static auto desc = MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + + int nmax; + parthenon::par_reduce( + DEFAULT_LOOP_PATTERN, "MaxSize", DevExecSpace(), 0, md->NumBlocks() - 1, + KOKKOS_LAMBDA(const int b, int &lmax) { + lmax = std::max(lmax, vmesh.GetSize(b, gas::cons::density())); + }, + Kokkos::Max(nmax)); + + // nothing to do + if (nmax == 1) return TaskStatus::complete; + + // Special cases + if (nmax == 2) { + return CoupleTwoFluids(md, dt); + } + + // General case + return CoupleNFluids(md, nmax, dt); +} + //---------------------------------------------------------------------------------------- //! template instantiations typedef Coordinates G; @@ -183,4 +235,18 @@ template TaskStatus DragSource(MD *md, const Real tt, const Real template TaskStatus DragSource(MD *md, const Real tt, const Real dt); template TaskStatus DragSource(MD *md, const Real tt, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); +template TaskStatus CoupleTwoFluids(MD *md, const Real dt); + +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); +template TaskStatus CoupleNFluids(MD *md, const int nmax, const Real dt); + } // namespace Drag diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index 4358cdd8..f26b2632 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -22,6 +22,7 @@ #include "utils/artemis_utils.hpp" #include "utils/diffusion/diffusion_coeff.hpp" #include "utils/eos/eos.hpp" +#include "utils/units.hpp" using namespace parthenon::package::prelude; using ArtemisUtils::EOS; @@ -53,8 +54,14 @@ namespace Drag { */ // ... Coupling types -enum class Coupling { simple_dust, self, null }; -// ... Drag models +enum class Coupling { simple_dust, self, full, null }; + +// ... Interspecies drag models for the full Chapman-Cowling coupling path +// hard_sphere - calibrated hard-sphere collision integrals (first backend) +// lj - Lennard-Jones surrogate collision integrals +enum class GasDragModel { hard_sphere, lj, null }; + +// ... Drag models for the simple dust-gas coupling path enum class DragModel { constant, stokes, null }; //---------------------------------------------------------------------------------------- @@ -65,6 +72,8 @@ inline Coupling ChooseDrag(const std::string choice) { return Coupling::self; } else if (choice == "simple_dust") { return Coupling::simple_dust; + } else if (choice == "full") { + return Coupling::full; } else { PARTHENON_FAIL("Bad choice of drag type"); return Coupling::null; @@ -158,420 +167,113 @@ struct StoppingTimeParams { }; //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Drag::SelfDragSourceImpl -//! \brief Implementation for self drag -template -TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt, - const Diffusion::DiffCoeffParams &dp, const EOS &eos_d, - const SelfDragParams &gasp, const SelfDragParams &dustp) { - PARTHENON_INSTRUMENT - using TE = parthenon::TopologicalElement; - auto pm = md->GetParentPointer(); - auto &resolved_pkgs = pm->resolved_packages; - - // Dimensionality - const int ndim = pm->ndim; - const int multi_d = (ndim >= 2); - const int three_d = (ndim == 3); - - // Extract artemis parameters - auto &artemis_pkg = pm->packages.Get("artemis"); - const bool do_gas = artemis_pkg->template Param("do_gas"); - const bool do_dust = artemis_pkg->template Param("do_dust"); - - // Extract gas parameters - Real de_switch = Null(); - Real dflr_gas = Null(); - Real sieflr_gas = Null(); - if (do_gas) { - auto &gas_pkg = pm->packages.Get("gas"); - de_switch = gas_pkg->template Param("de_switch"); - dflr_gas = gas_pkg->template Param("dfloor"); - sieflr_gas = gas_pkg->template Param("siefloor"); - } +//! \struct FullCouplingParams +//! \brief Parameters for the full Chapman-Cowling interspecies coupling +//! +//! Species are indexed 0..nspecies-1, matching gas.nspecies ordering. +//! mu_s[n] - molecular mass of species n in code-mass units +//! sigma_s[n] - effective hard-sphere diameter (collision cross-section radius) in +//! code +//! length units, or the LJ sigma parameter for the lj model +//! eps_s[n] - LJ epsilon/kB in EOS temperature units for the lj model +//! dof_s[n] - internal degrees of freedom per molecule (e.g. 3 for monatomic, +//! 5 for diatomic); used for energy-exchange weighting +struct FullCouplingParams { + GasDragModel model; + int nspecies; + Real kb_code; + ParArray1D mu_s; // molecular mass per species [code mass] + ParArray1D sigma_s; // collision diameter / LJ sigma per species [code length] + ParArray1D eps_s; // LJ epsilon/kB per species [EOS temperature units] + ParArray1D dof_s; // internal DOF per molecule per species + + FullCouplingParams() : model(GasDragModel::null), nspecies(0), kb_code(1.0) {} + + FullCouplingParams(ParameterInput *pin, const ArtemisUtils::Constants &constants, + ParArray1D gas_mu) { + const std::string block = "drag/full"; + nspecies = pin->GetOrAddInteger("gas", "nspecies", 1); + kb_code = constants.GetKBCode(); + + const std::string model_str = + pin->GetOrAddString(block, "collision_model", "hard_sphere"); + if (model_str == "hard_sphere") { + model = GasDragModel::hard_sphere; + } else if (model_str == "lj") { + model = GasDragModel::lj; + } else { + PARTHENON_FAIL("drag/full collision_model must be hard_sphere or lj"); + model = GasDragModel::null; + } - // Extract dust parameters - Real dflr_dust = Null(); - if (do_dust) { - auto &dust_pkg = pm->packages.Get("dust"); - dflr_dust = dust_pkg->template Param("dfloor"); - } - const auto &cpars = artemis_pkg->template Param("coord_params"); - - const Real x1min = artemis_pkg->template Param("x1min"); - const Real x1max = artemis_pkg->template Param("x1max"); - const Real x2min = artemis_pkg->template Param("x2min"); - const Real x2max = artemis_pkg->template Param("x2max"); - const Real x3min = artemis_pkg->template Param("x3min"); - const Real x3max = artemis_pkg->template Param("x3max"); - - // Packing and indexing - static auto desc = - MakePackDescriptor(resolved_pkgs.get()); - auto vmesh = desc.GetPack(md); - static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); - auto vg = desc_g.GetPack(md); - const auto ib = md->GetBoundsI(IndexDomain::interior); - const auto jb = md->GetBoundsJ(IndexDomain::interior); - const auto kb = md->GetBoundsK(IndexDomain::interior); - - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "SelfDrag", parthenon::DevExecSpace(), 0, md->NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { - // Extract coordinates - geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); - const auto &xv = coords.GetCellCenter(vg, b, k, j, i); - const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); - const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); - - // Compute the (gas) ramp for this cell - // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) - if (do_gas) { - const Real fx1 = - dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * - SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + - gasp.orate[0] * ((xv[0] > gasp.ox[0]) * - SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max)))); - const Real fx2 = - multi_d * dt * - (gasp.irate[1] * ((xv[1] < gasp.ix[1]) * - SQR((xv[1] - gasp.ix[1]) / (gasp.ix[1] - x2min))) + - gasp.orate[1] * ((xv[1] > gasp.ox[1]) * - SQR((xv[1] - gasp.ox[1]) / (gasp.ox[1] - x2max)))); - const Real fx3 = - three_d * dt * - (gasp.irate[2] * ((xv[2] < gasp.ix[2]) * - SQR((xv[2] - gasp.ix[2]) / (gasp.ix[2] - x3min))) + - gasp.orate[2] * ((xv[2] > gasp.ox[2]) * - SQR((xv[2] - gasp.ox[2]) / (gasp.ox[2] - x3max)))); - - // Update gas momenta and total energy - for (int n = 0; n < vmesh.GetSize(b, gas::cons::density()); ++n) { - // Extract state vector - Real &dens = vmesh(b, gas::cons::density(n), k, j, i); - Real &mom1 = vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i); - Real &mom2 = vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i); - Real &mom3 = vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i); - Real &etot = vmesh(b, gas::cons::total_energy(n), k, j, i); - - // Apply density floor - const bool dfloor = (dens > dflr_gas); - dens = (dfloor)*dens + (!dfloor) * dflr_gas; - - // Compute SIE via dual energy formalism and apply floor - Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); - const Real efloor = (sieg > sieflr_gas); - sieg = (efloor)*sieg + (!efloor) * sieflr_gas; - - // Get diffusion coefficient - Diffusion::DiffusionCoeff dcoeff; - const Real mu = dcoeff.Get(dp, coords, xv, dens, sieg, eos_d); - const Real vR = -1.5 * mu / (xcyl[0] * dens); - const Real vg[3] = {mom1 / (hx[0] * dens), mom2 / (hx[1] * dens), - mom3 / (hx[2] * dens)}; - const Real vd[3] = {ex1[0] * vR, ex2[0] * vR, ex3[0] * vR}; - - // Apply update - // Ep - E = 0.5 d ( vp^2 - v^2 ) - // (vp-v) . (vp + v) = dv . (2v + dv) = 2 dv.v + dv.dv - const Real dm1 = -fx1 * dens * (vg[0] - vd[0]) / (1.0 + fx1); - const Real dm2 = -fx2 * dens * (vg[1] - vd[1]) / (1.0 + fx2); - const Real dm3 = -fx3 * dens * (vg[2] - vd[2]) / (1.0 + fx3); - mom1 += hx[0] * dm1; - mom2 += hx[1] * dm2; - mom3 += hx[2] * dm3; - etot += dm1 * (vg[0] + 0.5 * dm1 / dens) + dm2 * (vg[1] + 0.5 * dm2 / dens) + - dm3 * (vg[2] + 0.5 * dm3 / dens); - - // Apply total energy floor - const Real utmp = dens * sieflr_gas; - etot = std::max(etot, utmp); - } - } - - // Compute the (dust) ramp for this cell - if (do_dust) { - const Real fx1 = - dt * - (dustp.irate[0] * ((xv[0] < dustp.ix[0]) * - SQR((xv[0] - dustp.ix[0]) / (dustp.ix[0] - x1min))) + - dustp.orate[0] * ((xv[0] > dustp.ox[0]) * - SQR((xv[0] - dustp.ox[0]) / (dustp.ox[0] - x1max)))); - const Real fx2 = - multi_d * dt * - (dustp.irate[1] * ((xv[1] < dustp.ix[1]) * - SQR((xv[1] - dustp.ix[1]) / (dustp.ix[1] - x2min))) + - dustp.orate[1] * ((xv[1] > dustp.ox[1]) * - SQR((xv[1] - dustp.ox[1]) / (dustp.ox[1] - x2max)))); - const Real fx3 = - three_d * dt * - (dustp.irate[2] * ((xv[2] < dustp.ix[2]) * - SQR((xv[2] - dustp.ix[2]) / (dustp.ix[2] - x3min))) + - dustp.orate[2] * ((xv[2] > dustp.ox[2]) * - SQR((xv[2] - dustp.ox[2]) / (dustp.ox[2] - x3max)))); - - // Update dust momenta - for (int n = 0; n < vmesh.GetSize(b, dust::cons::density()); ++n) { - // Extract state vector - Real &dens = vmesh(b, dust::cons::density(n), k, j, i); - Real &mom1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i); - Real &mom2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i); - Real &mom3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i); - - // Apply density floor - const bool dfloor = (dens > dflr_dust); - dens = (dfloor)*dens + (!dfloor) * dflr_dust; - - // Apply update - mom1 -= fx1 * mom1 / (1.0 + fx1); - mom2 -= fx2 * mom2 / (1.0 + fx2); - mom3 -= fx3 * mom3 / (1.0 + fx3); - } - } - }); - - return TaskStatus::complete; -} + mu_s = ParArray1D("drag_mu", nspecies); + sigma_s = ParArray1D("drag_sigma", nspecies); + eps_s = ParArray1D("drag_eps", nspecies); + dof_s = ParArray1D("drag_dof", nspecies); + + auto h_mu = mu_s.GetHostMirror(); + auto h_sigma = sigma_s.GetHostMirror(); + auto h_eps = eps_s.GetHostMirror(); + auto h_dof = dof_s.GetHostMirror(); + const auto gas_mu_h = gas_mu.GetHostMirrorAndCopy(); + + std::vector mu_v = pin->GetVector(block, "mu"); + std::vector sigma_v = pin->GetVector(block, "sigma"); + + // Optional: LJ epsilon and internal DOF (default to monatomic hard-sphere) + std::vector eps_v(nspecies, 0.0); + std::vector dof_v(nspecies, 3.0); + const bool has_eps = pin->DoesParameterExist(block, "eps"); + if (has_eps) eps_v = pin->GetVector(block, "eps"); + if (pin->DoesParameterExist(block, "dof")) dof_v = pin->GetVector(block, "dof"); + + PARTHENON_REQUIRE(static_cast(mu_v.size()) == nspecies, + "drag/full mu must have nspecies entries"); + PARTHENON_REQUIRE(static_cast(sigma_v.size()) == nspecies, + "drag/full sigma must have nspecies entries"); + PARTHENON_REQUIRE(static_cast(dof_v.size()) == nspecies, + "drag/full dof must have nspecies entries"); + if (model == GasDragModel::lj) { + PARTHENON_REQUIRE(has_eps, "drag/full eps is required for collision_model = lj"); + PARTHENON_REQUIRE(static_cast(eps_v.size()) == nspecies, + "drag/full eps must have nspecies entries"); + } else { + PARTHENON_REQUIRE(!has_eps, "drag/full eps is only valid for collision_model = lj"); + } -//---------------------------------------------------------------------------------------- -//! \fn TaskStatus Drag::SimpleDragSourceImpl -//! \brief Implementation for simple drag -template -TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real dt, - const Diffusion::DiffCoeffParams &dp, const EOS &eos_d, - const SelfDragParams &gasp, const SelfDragParams &dustp, - const StoppingTimeParams &tp) { - PARTHENON_INSTRUMENT - using TE = parthenon::TopologicalElement; - auto pm = md->GetParentPointer(); - auto &resolved_pkgs = pm->resolved_packages; - - // Dimensionality - const int ndim = pm->ndim; - const int multi_d = (ndim >= 2); - const int three_d = (ndim == 3); - auto &artemis_pkg = pm->packages.Get("artemis"); - - // Extract gas package and params - auto &gas_pkg = pm->packages.Get("gas"); - const Real de_switch = gas_pkg->template Param("de_switch"); - const Real dflr_gas = gas_pkg->template Param("dfloor"); - const Real sieflr_gas = gas_pkg->template Param("siefloor"); - - const Real x1min = artemis_pkg->template Param("x1min"); - const Real x1max = artemis_pkg->template Param("x1max"); - const Real x2min = artemis_pkg->template Param("x2min"); - const Real x2max = artemis_pkg->template Param("x2max"); - const Real x3min = artemis_pkg->template Param("x3min"); - const Real x3max = artemis_pkg->template Param("x3max"); - - // Extract dust package and params - auto &dust_pkg = pm->packages.Get("dust"); - const auto &sizes = dust_pkg->template Param>("sizes"); - const auto grain_density = dust_pkg->template Param("grain_density"); - const Real dflr_dust = dust_pkg->template Param("dfloor"); - - const auto &cpars = artemis_pkg->template Param("coord_params"); - - // Packing and indexing - static auto desc = - MakePackDescriptor(resolved_pkgs.get()); - auto vmesh = desc.GetPack(md); - static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); - auto vg = desc_g.GetPack(md); - const auto ib = md->GetBoundsI(IndexDomain::interior); - const auto jb = md->GetBoundsJ(IndexDomain::interior); - const auto kb = md->GetBoundsK(IndexDomain::interior); - - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "SimpleDrag", parthenon::DevExecSpace(), 0, - md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { - // Extract coordinates - geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); - const auto &xv = coords.GetCellCenter(vg, b, k, j, i); - const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); - const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); - - // Compute the ramp for this cell - // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) - const std::array bg{ - dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * - SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + - gasp.orate[0] * ((xv[0] > gasp.ox[0]) * - SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max)))), - multi_d * dt * - (gasp.irate[1] * ((xv[1] < gasp.ix[1]) * - SQR((xv[1] - gasp.ix[1]) / (gasp.ix[1] - x2min))) + - gasp.orate[1] * ((xv[1] > gasp.ox[1]) * - SQR((xv[1] - gasp.ox[1]) / (gasp.ox[1] - x2max)))), - three_d * dt * - (gasp.irate[2] * ((xv[2] < gasp.ix[2]) * - SQR((xv[2] - gasp.ix[2]) / (gasp.ix[2] - x3min))) + - gasp.orate[2] * ((xv[2] > gasp.ox[2]) * - SQR((xv[2] - gasp.ox[2]) / (gasp.ox[2] - x3max))))}; - const std::array bd{ - dt * (dustp.irate[0] * ((xv[0] < dustp.ix[0]) * - SQR((xv[0] - dustp.ix[0]) / (dustp.ix[0] - x1min))) + - dustp.orate[0] * ((xv[0] > dustp.ox[0]) * - SQR((xv[0] - dustp.ox[0]) / (dustp.ox[0] - x1max)))), - multi_d * dt * - (dustp.irate[1] * ((xv[1] < dustp.ix[1]) * - SQR((xv[1] - dustp.ix[1]) / (dustp.ix[1] - x2min))) + - dustp.orate[1] * ((xv[1] > dustp.ox[1]) * - SQR((xv[1] - dustp.ox[1]) / (dustp.ox[1] - x2max)))), - three_d * dt * - (dustp.irate[2] * ((xv[2] < dustp.ix[2]) * - SQR((xv[2] - dustp.ix[2]) / (dustp.ix[2] - x3min))) + - dustp.orate[2] * ((xv[2] > dustp.ox[2]) * - SQR((xv[2] - dustp.ox[2]) / (dustp.ox[2] - x3max))))}; - - // Extract gas state vector - // NOTE(@pdmullen): Assumes single gas species - Real &dg = vmesh(b, gas::cons::density(0), k, j, i); - Real &gmom1 = vmesh(b, gas::cons::momentum(VI(0, 0)), k, j, i); - Real &gmom2 = vmesh(b, gas::cons::momentum(VI(0, 1)), k, j, i); - Real &gmom3 = vmesh(b, gas::cons::momentum(VI(0, 2)), k, j, i); - Real &etot = vmesh(b, gas::cons::total_energy(0), k, j, i); - - // Apply density floor - const bool dfloor = (dg > dflr_gas); - dg = (dfloor)*dg + (!dfloor) * dflr_gas; - - // Compute SIE via dual energy formalism and apply floor - Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, 0, k, j, i, de_switch, hx); - const Real efloor = (sieg > sieflr_gas); - sieg = (efloor)*sieg + (!efloor) * sieflr_gas; - - // Stash gas velocity - const std::array vg{gmom1 / (hx[0] * dg), gmom2 / (hx[1] * dg), - gmom3 / (hx[2] * dg)}; - - // Target gas velocity - Diffusion::DiffusionCoeff dcoeff; - const Real mu = dcoeff.Get(dp, coords, xv, dg, sieg, eos_d); - const Real vR = -1.5 * mu / (xcyl[0] * dg); - const std::array vt{ex1[0] * vR, ex2[0] * vR, ex3[0] * vR}; - - // Extract Stokes specific parameters - [[maybe_unused]] auto &grain_density_ = grain_density; - [[maybe_unused]] Real vth = Null(); - if constexpr (DRAG == DragModel::stokes) { - const Real gm1 = eos_d.GruneisenParamFromDensityInternalEnergy(dg, sieg); - vth = std::sqrt(8.0 / M_PI * gm1 * sieg); - } - - // First pass to collect \sum rho' and \sum rho' v and compute new vg - std::array fd{0.0, 0.0, 0.0}; - std::array fvd{0.0, 0.0, 0.0}; - std::array vdt{0.0, 0.0, 0.0}; - const int nspecies = vmesh.GetSize(b, dust::cons::density()); - for (int n = 0; n < nspecies; ++n) { - // Extract state vector - Real &dens = vmesh(b, dust::cons::density(n), k, j, i); - Real &dmom1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i); - Real &dmom2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i); - Real &dmom3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i); - - // Apply density floor - const bool dfloor = (dens > dflr_dust); - dens = (dfloor)*dens + (!dfloor) * dflr_dust; - - // Stash nth dust velocity - const std::array vd{dmom1 / (hx[0] * dens), dmom2 / (hx[1] * dens), - dmom3 / (hx[2] * dens)}; - - // Coupling - const auto id = vmesh(b, dust::cons::density(n)).sparse_id; - Real tc = tp.tau(id); - [[maybe_unused]] auto &sizes_ = sizes; - if constexpr (DRAG == DragModel::stokes) { - tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / - dg * sizes_(id) / vth)); - } - const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); - for (int d = 0; d < 3; d++) { - const Real rhop = dens * alpha / (1.0 + alpha + bd[d]); - fd[d] += rhop * (1.0 + bd[d]); - fvd[d] += rhop * (vd[d] + bd[d] * vdt[d]); - } - } - - // New vgas - std::array vgp{Null(), Null(), Null()}; - for (int d = 0; d < 3; d++) { - vgp[d] = (dg * (vg[d] + bg[d] * vt[d]) + fvd[d]) / (dg * (1.0 + bg[d]) + fd[d]); - } - - // Second pass to update all momenta - fvd = {0.0, 0.0, 0.0}; - std::array delta_g{0.0, 0.0, 0.0}; - for (int n = 0; n < nspecies; ++n) { - // Extract dust density - // NOTE(@pdmullen): Dust density already floored above - const Real &dens = vmesh(b, dust::cons::density(n), k, j, i); - - // Stash nth dust velocity - const std::array vd{ - vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * dens), - vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * dens), - vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * dens)}; - - // Coupling - const auto id = vmesh(b, dust::cons::density(n)).sparse_id; - Real tc = tp.tau(id); - [[maybe_unused]] auto &sizes_ = sizes; - if constexpr (DRAG == DragModel::stokes) { - tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / - dg * sizes_(id) / vth)); - } - const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); - // Update dust momenta - for (int d = 0; d < 3; d++) { - Real delta_d = 0.; - const Real rhop = dens * alpha / (1.0 + alpha + bd[d]); - const Real delta = rhop * ((vgp[d] - vd[d] + bd[d] * (vgp[d] - vdt[d]))); - delta_d += delta; - delta_g[d] -= delta; - - // self-drag coupling - delta_d -= bd[d] * dens / (1. + alpha + bd[d]) * - (vd[d] - vdt[d] + alpha * (vgp[d] - vdt[d])); - fvd[d] += rhop * (vd[d] - vt[d] + bd[d] * (vdt[d] - vt[d])); - vmesh(b, dust::cons::momentum(VI(n, d)), k, j, i) += hx[d] * delta_d; - } - } - - // Final update to gas momenta and total energy - for (int d = 0; d < 3; d++) { - const Real prefac = dg * bg[d] / (1.0 + bg[d] + fd[d]); - delta_g[d] -= prefac * (dg * (vg[d] - vt[d]) + fvd[d]); - const Real vn = vg[d] + delta_g[d] / dg; - vmesh(b, gas::cons::momentum(VI(0, d)), k, j, i) += hx[d] * delta_g[d]; - etot += 0.5 * (vg[d] + vn) * delta_g[d]; - } - - // Apply total energy floors - const Real utmp = dg * sieflr_gas; - etot = std::max(etot, utmp); - }); - - return TaskStatus::complete; -} + for (int n = 0; n < nspecies; ++n) { + PARTHENON_REQUIRE(mu_v[n] > 0.0, "drag/full mu entries must be positive"); + PARTHENON_REQUIRE(sigma_v[n] > 0.0, "drag/full sigma entries must be positive"); + PARTHENON_REQUIRE(dof_v[n] > 0.0, "drag/full dof entries must be positive"); + PARTHENON_REQUIRE(std::abs(mu_v[n] - gas_mu_h(n)) <= + 1.0e-12 * std::abs(gas_mu_h(n)), + "drag/full mu must match gas mu"); + if (model == GasDragModel::lj) { + PARTHENON_REQUIRE(eps_v[n] > 0.0, "drag/full eps entries must be positive"); + } + h_mu(n) = mu_v[n] * constants.GetAMUCode(); + h_sigma(n) = sigma_v[n]; + h_eps(n) = eps_v[n]; + h_dof(n) = dof_v[n]; + } + mu_s.DeepCopy(h_mu); + sigma_s.DeepCopy(h_sigma); + eps_s.DeepCopy(h_eps); + dof_s.DeepCopy(h_dof); + } +}; //---------------------------------------------------------------------------------------- -std::shared_ptr Initialize(ParameterInput *pin); +std::shared_ptr Initialize(ParameterInput *pin, + const ArtemisUtils::Constants &constants, + Packages_t &packages); template TaskStatus DragSource(MeshData *md, const Real time, const Real dt); +template +TaskStatus ApplyClosure(MeshData *md, const Real dt); + } // namespace Drag #endif // DRAG_DRAG_HPP_ diff --git a/src/drag/drag_impl.hpp b/src/drag/drag_impl.hpp new file mode 100644 index 00000000..7d0ec60a --- /dev/null +++ b/src/drag/drag_impl.hpp @@ -0,0 +1,884 @@ +//======================================================================================== +// (C) (or copyright) 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// icense in this material to reproduce, prepare derivative works, distribute copies to +// the pubic, perform pubicly and display pubicly, and to permit others to do so. +//======================================================================================== +#ifndef DRAG_DRAG_IMPL_HPP_ +#define DRAG_DRAG_IMPL_HPP_ + +// Parthenon includes +#include + +// KokkosKernels batched LU +#include +#include + +// Artemis includes +#include "artemis.hpp" +#include "collision_integrals.hpp" +#include "drag.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/diffusion/diffusion_coeff.hpp" +#include "utils/eos/eos.hpp" + +using namespace parthenon::package::prelude; +using ArtemisUtils::EOS; +using ArtemisUtils::VI; + +namespace Drag { + +KOKKOS_FORCEINLINE_FUNCTION Real ClampRoundoffHeat(const Real q, const Real scale) { + const Real tol = 1.0e-10 * (std::abs(scale) + 1.0); + PARTHENON_DEBUG_REQUIRE(q >= -tol, "Gas-gas coupling generated negative drift heating"); + return (q < 0.0 && std::abs(q) <= tol) ? 0.0 : q; +} + +KOKKOS_FORCEINLINE_FUNCTION void ApplyEnergyFloor(const Real rho, const Real sie_floor, + const Real ke, Real &etot, Real &eint) { + eint = std::max(eint, rho * sie_floor); + etot = std::max(etot, eint + ke); +} + +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus Drag::SelfDragSourceImpl +//! \brief Implementation for self drag +template +TaskStatus SelfDragSourceImpl(MeshData *md, const Real time, const Real dt, + const Diffusion::DiffCoeffParams &dp, const EOS &eos_d, + const SelfDragParams &gasp, const SelfDragParams &dustp) { + PARTHENON_INSTRUMENT + using TE = parthenon::TopologicalElement; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Dimensionaity + const int ndim = pm->ndim; + const int multi_d = (ndim >= 2); + const int three_d = (ndim == 3); + + // Extract artemis parameters + auto &artemis_pkg = pm->packages.Get("artemis"); + const bool do_gas = artemis_pkg->template Param("do_gas"); + const bool do_dust = artemis_pkg->template Param("do_dust"); + + // Extract gas parameters + Real de_switch = Null(); + Real dflr_gas = Null(); + Real sieflr_gas = Null(); + if (do_gas) { + auto &gas_pkg = pm->packages.Get("gas"); + de_switch = gas_pkg->template Param("de_switch"); + dflr_gas = gas_pkg->template Param("dfloor"); + sieflr_gas = gas_pkg->template Param("siefloor"); + } + + // Extract dust parameters + Real dflr_dust = Null(); + if (do_dust) { + auto &dust_pkg = pm->packages.Get("dust"); + dflr_dust = dust_pkg->template Param("dfloor"); + } + const auto &cpars = artemis_pkg->template Param("coord_params"); + + const Real x1min = artemis_pkg->template Param("x1min"); + const Real x1max = artemis_pkg->template Param("x1max"); + const Real x2min = artemis_pkg->template Param("x2min"); + const Real x2max = artemis_pkg->template Param("x2max"); + const Real x3min = artemis_pkg->template Param("x3min"); + const Real x3max = artemis_pkg->template Param("x3max"); + + // Packing and indexing + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SelfDrag", parthenon::DevExecSpace(), 0, md->NumBlocks() - 1, + kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + // Extract coordinates + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &xv = coords.GetCellCenter(vg, b, k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + + // Compute the (gas) ramp for this cell + // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) + if (do_gas) { + const Real fx1 = + dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * + SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + + gasp.orate[0] * ((xv[0] > gasp.ox[0]) * + SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max)))); + const Real fx2 = + multi_d * dt * + (gasp.irate[1] * ((xv[1] < gasp.ix[1]) * + SQR((xv[1] - gasp.ix[1]) / (gasp.ix[1] - x2min))) + + gasp.orate[1] * ((xv[1] > gasp.ox[1]) * + SQR((xv[1] - gasp.ox[1]) / (gasp.ox[1] - x2max)))); + const Real fx3 = + three_d * dt * + (gasp.irate[2] * ((xv[2] < gasp.ix[2]) * + SQR((xv[2] - gasp.ix[2]) / (gasp.ix[2] - x3min))) + + gasp.orate[2] * ((xv[2] > gasp.ox[2]) * + SQR((xv[2] - gasp.ox[2]) / (gasp.ox[2] - x3max)))); + + // Update gas momenta and total energy + for (int n = 0; n < vmesh.GetSize(b, gas::cons::density()); ++n) { + // Extract state vector + Real &dens = vmesh(b, gas::cons::density(n), k, j, i); + Real &mom1 = vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i); + Real &mom2 = vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i); + Real &mom3 = vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i); + Real &etot = vmesh(b, gas::cons::total_energy(n), k, j, i); + + // Apply density floor + const bool dfloor = (dens > dflr_gas); + dens = (dfloor)*dens + (!dfloor) * dflr_gas; + + // Compute SIE via dual energy formaism and apply floor + Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + const Real efloor = (sieg > sieflr_gas); + sieg = (efloor)*sieg + (!efloor) * sieflr_gas; + + // Get diffusion coefficient + Diffusion::DiffusionCoeff dcoeff; + const Real mu = dcoeff.Get(dp, coords, xv, dens, sieg, eos_d); + const Real vR = -1.5 * mu / (xcyl[0] * dens); + const Real vg[3] = {mom1 / (hx[0] * dens), mom2 / (hx[1] * dens), + mom3 / (hx[2] * dens)}; + const Real vd[3] = {ex1[0] * vR, ex2[0] * vR, ex3[0] * vR}; + + // Apply update + // Ep - E = 0.5 d ( vp^2 - v^2 ) + // (vp-v) . (vp + v) = dv . (2v + dv) = 2 dv.v + dv.dv + const Real dm1 = -fx1 * dens * (vg[0] - vd[0]) / (1.0 + fx1); + const Real dm2 = -fx2 * dens * (vg[1] - vd[1]) / (1.0 + fx2); + const Real dm3 = -fx3 * dens * (vg[2] - vd[2]) / (1.0 + fx3); + mom1 += hx[0] * dm1; + mom2 += hx[1] * dm2; + mom3 += hx[2] * dm3; + etot += dm1 * (vg[0] + 0.5 * dm1 / dens) + dm2 * (vg[1] + 0.5 * dm2 / dens) + + dm3 * (vg[2] + 0.5 * dm3 / dens); + + // Apply total energy floor + const Real utmp = dens * sieflr_gas; + etot = std::max(etot, utmp); + } + } + + // Compute the (dust) ramp for this cell + if (do_dust) { + const Real fx1 = + dt * + (dustp.irate[0] * ((xv[0] < dustp.ix[0]) * + SQR((xv[0] - dustp.ix[0]) / (dustp.ix[0] - x1min))) + + dustp.orate[0] * ((xv[0] > dustp.ox[0]) * + SQR((xv[0] - dustp.ox[0]) / (dustp.ox[0] - x1max)))); + const Real fx2 = + multi_d * dt * + (dustp.irate[1] * ((xv[1] < dustp.ix[1]) * + SQR((xv[1] - dustp.ix[1]) / (dustp.ix[1] - x2min))) + + dustp.orate[1] * ((xv[1] > dustp.ox[1]) * + SQR((xv[1] - dustp.ox[1]) / (dustp.ox[1] - x2max)))); + const Real fx3 = + three_d * dt * + (dustp.irate[2] * ((xv[2] < dustp.ix[2]) * + SQR((xv[2] - dustp.ix[2]) / (dustp.ix[2] - x3min))) + + dustp.orate[2] * ((xv[2] > dustp.ox[2]) * + SQR((xv[2] - dustp.ox[2]) / (dustp.ox[2] - x3max)))); + + // Update dust momenta + for (int n = 0; n < vmesh.GetSize(b, dust::cons::density()); ++n) { + // Extract state vector + Real &dens = vmesh(b, dust::cons::density(n), k, j, i); + Real &mom1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i); + Real &mom2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i); + Real &mom3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i); + + // Apply density floor + const bool dfloor = (dens > dflr_dust); + dens = (dfloor)*dens + (!dfloor) * dflr_dust; + + // Apply update + mom1 -= fx1 * mom1 / (1.0 + fx1); + mom2 -= fx2 * mom2 / (1.0 + fx2); + mom3 -= fx3 * mom3 / (1.0 + fx3); + } + } + }); + + return TaskStatus::complete; +} + +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus Drag::SimpleDragSourceImpl +//! \brief Implementation for simple drag +template +TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real dt, + const Diffusion::DiffCoeffParams &dp, const EOS &eos_d, + const SelfDragParams &gasp, const SelfDragParams &dustp, + const StoppingTimeParams &tp) { + PARTHENON_INSTRUMENT + using TE = parthenon::TopologicalElement; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Dimensionaity + const int ndim = pm->ndim; + const int multi_d = (ndim >= 2); + const int three_d = (ndim == 3); + auto &artemis_pkg = pm->packages.Get("artemis"); + + // Extract gas package and params + auto &gas_pkg = pm->packages.Get("gas"); + const Real de_switch = gas_pkg->template Param("de_switch"); + const Real dflr_gas = gas_pkg->template Param("dfloor"); + const Real sieflr_gas = gas_pkg->template Param("siefloor"); + + const Real x1min = artemis_pkg->template Param("x1min"); + const Real x1max = artemis_pkg->template Param("x1max"); + const Real x2min = artemis_pkg->template Param("x2min"); + const Real x2max = artemis_pkg->template Param("x2max"); + const Real x3min = artemis_pkg->template Param("x3min"); + const Real x3max = artemis_pkg->template Param("x3max"); + + // Extract dust package and params + auto &dust_pkg = pm->packages.Get("dust"); + const auto &sizes = dust_pkg->template Param>("sizes"); + const auto grain_density = dust_pkg->template Param("grain_density"); + const Real dflr_dust = dust_pkg->template Param("dfloor"); + + const auto &cpars = artemis_pkg->template Param("coord_params"); + + // Packing and indexing + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SimpleDrag", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + // Extract coordinates + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &xv = coords.GetCellCenter(vg, b, k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + + // Compute the ramp for this cell + // Ramps are quadratic, eg. the left regions is SQR( (X - ix)/(ix - xmin) ) + const std::array bg{ + dt * (gasp.irate[0] * ((xv[0] < gasp.ix[0]) * + SQR((xv[0] - gasp.ix[0]) / (gasp.ix[0] - x1min))) + + gasp.orate[0] * ((xv[0] > gasp.ox[0]) * + SQR((xv[0] - gasp.ox[0]) / (gasp.ox[0] - x1max)))), + multi_d * dt * + (gasp.irate[1] * ((xv[1] < gasp.ix[1]) * + SQR((xv[1] - gasp.ix[1]) / (gasp.ix[1] - x2min))) + + gasp.orate[1] * ((xv[1] > gasp.ox[1]) * + SQR((xv[1] - gasp.ox[1]) / (gasp.ox[1] - x2max)))), + three_d * dt * + (gasp.irate[2] * ((xv[2] < gasp.ix[2]) * + SQR((xv[2] - gasp.ix[2]) / (gasp.ix[2] - x3min))) + + gasp.orate[2] * ((xv[2] > gasp.ox[2]) * + SQR((xv[2] - gasp.ox[2]) / (gasp.ox[2] - x3max))))}; + const std::array bd{ + dt * (dustp.irate[0] * ((xv[0] < dustp.ix[0]) * + SQR((xv[0] - dustp.ix[0]) / (dustp.ix[0] - x1min))) + + dustp.orate[0] * ((xv[0] > dustp.ox[0]) * + SQR((xv[0] - dustp.ox[0]) / (dustp.ox[0] - x1max)))), + multi_d * dt * + (dustp.irate[1] * ((xv[1] < dustp.ix[1]) * + SQR((xv[1] - dustp.ix[1]) / (dustp.ix[1] - x2min))) + + dustp.orate[1] * ((xv[1] > dustp.ox[1]) * + SQR((xv[1] - dustp.ox[1]) / (dustp.ox[1] - x2max)))), + three_d * dt * + (dustp.irate[2] * ((xv[2] < dustp.ix[2]) * + SQR((xv[2] - dustp.ix[2]) / (dustp.ix[2] - x3min))) + + dustp.orate[2] * ((xv[2] > dustp.ox[2]) * + SQR((xv[2] - dustp.ox[2]) / (dustp.ox[2] - x3max))))}; + + // Extract gas state vector + // NOTE(@pdmullen): Assumes single gas species + Real &dg = vmesh(b, gas::cons::density(0), k, j, i); + Real &gmom1 = vmesh(b, gas::cons::momentum(VI(0, 0)), k, j, i); + Real &gmom2 = vmesh(b, gas::cons::momentum(VI(0, 1)), k, j, i); + Real &gmom3 = vmesh(b, gas::cons::momentum(VI(0, 2)), k, j, i); + Real &etot = vmesh(b, gas::cons::total_energy(0), k, j, i); + + // Apply density floor + const bool dfloor = (dg > dflr_gas); + dg = (dfloor)*dg + (!dfloor) * dflr_gas; + + // Compute SIE via dual energy formaism and apply floor + Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, 0, k, j, i, de_switch, hx); + const Real efloor = (sieg > sieflr_gas); + sieg = (efloor)*sieg + (!efloor) * sieflr_gas; + + // Stash gas velocity + const std::array vg{gmom1 / (hx[0] * dg), gmom2 / (hx[1] * dg), + gmom3 / (hx[2] * dg)}; + + // Target gas velocity + Diffusion::DiffusionCoeff dcoeff; + const Real mu = dcoeff.Get(dp, coords, xv, dg, sieg, eos_d); + const Real vR = -1.5 * mu / (xcyl[0] * dg); + const std::array vt{ex1[0] * vR, ex2[0] * vR, ex3[0] * vR}; + + // Extract Stokes specific parameters + [[maybe_unused]] auto &grain_density_ = grain_density; + [[maybe_unused]] Real vth = Null(); + if constexpr (DRAG == DragModel::stokes) { + const Real gm1 = eos_d.GruneisenParamFromDensityInternalEnergy(dg, sieg); + vth = std::sqrt(8.0 / M_PI * gm1 * sieg); + } + + // First pass to collect \sum rho' and \sum rho' v and compute new vg + std::array fd{0.0, 0.0, 0.0}; + std::array fvd{0.0, 0.0, 0.0}; + std::array vdt{0.0, 0.0, 0.0}; + const int nspecies = vmesh.GetSize(b, dust::cons::density()); + for (int n = 0; n < nspecies; ++n) { + // Extract state vector + Real &dens = vmesh(b, dust::cons::density(n), k, j, i); + Real &dmom1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i); + Real &dmom2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i); + Real &dmom3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i); + + // Apply density floor + const bool dfloor = (dens > dflr_dust); + dens = (dfloor)*dens + (!dfloor) * dflr_dust; + + // Stash nth dust velocity + const std::array vd{dmom1 / (hx[0] * dens), dmom2 / (hx[1] * dens), + dmom3 / (hx[2] * dens)}; + + // Coupling + const auto id = vmesh(b, dust::cons::density(n)).sparse_id; + Real tc = tp.tau(id); + [[maybe_unused]] auto &sizes_ = sizes; + if constexpr (DRAG == DragModel::stokes) { + tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / + dg * sizes_(id) / vth)); + } + const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); + for (int d = 0; d < 3; d++) { + const Real rhop = dens * alpha / (1.0 + alpha + bd[d]); + fd[d] += rhop * (1.0 + bd[d]); + fvd[d] += rhop * (vd[d] + bd[d] * vdt[d]); + } + } + + // New vgas + std::array vgp{Null(), Null(), Null()}; + for (int d = 0; d < 3; d++) { + vgp[d] = (dg * (vg[d] + bg[d] * vt[d]) + fvd[d]) / (dg * (1.0 + bg[d]) + fd[d]); + } + + // Second pass to update all momenta + fvd = {0.0, 0.0, 0.0}; + std::array delta_g{0.0, 0.0, 0.0}; + for (int n = 0; n < nspecies; ++n) { + // Extract dust density + // NOTE(@pdmullen): Dust density already floored above + const Real &dens = vmesh(b, dust::cons::density(n), k, j, i); + + // Stash nth dust velocity + const std::array vd{ + vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * dens), + vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * dens), + vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * dens)}; + + // Coupling + const auto id = vmesh(b, dust::cons::density(n)).sparse_id; + Real tc = tp.tau(id); + [[maybe_unused]] auto &sizes_ = sizes; + if constexpr (DRAG == DragModel::stokes) { + tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / + dg * sizes_(id) / vth)); + } + const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); + // Update dust momenta + for (int d = 0; d < 3; d++) { + Real delta_d = 0.; + const Real rhop = dens * alpha / (1.0 + alpha + bd[d]); + const Real delta = rhop * ((vgp[d] - vd[d] + bd[d] * (vgp[d] - vdt[d]))); + delta_d += delta; + delta_g[d] -= delta; + + // self-drag coupling + delta_d -= bd[d] * dens / (1. + alpha + bd[d]) * + (vd[d] - vdt[d] + alpha * (vgp[d] - vdt[d])); + fvd[d] += rhop * (vd[d] - vt[d] + bd[d] * (vdt[d] - vt[d])); + vmesh(b, dust::cons::momentum(VI(n, d)), k, j, i) += hx[d] * delta_d; + } + } + + // Final update to gas momenta and total energy + for (int d = 0; d < 3; d++) { + const Real prefac = dg * bg[d] / (1.0 + bg[d] + fd[d]); + delta_g[d] -= prefac * (dg * (vg[d] - vt[d]) + fvd[d]); + const Real vn = vg[d] + delta_g[d] / dg; + vmesh(b, gas::cons::momentum(VI(0, d)), k, j, i) += hx[d] * delta_g[d]; + etot += 0.5 * (vg[d] + vn) * delta_g[d]; + } + + // Apply total energy floors + const Real utmp = dg * sieflr_gas; + etot = std::max(etot, utmp); + }); + + return TaskStatus::complete; +} + +template +TaskStatus CoupleTwoFluids(MeshData *md, const Real dt) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; + + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Extract gas package and params + auto &gas_pkg = pm->packages.Get("gas"); + const auto eos_d = gas_pkg->template Param>("eos_d"); + const auto dflr_gas = gas_pkg->template Param("dfloor"); + const auto sieflr_gas = gas_pkg->template Param("siefloor"); + const auto de_switch = gas_pkg->template Param("de_switch"); + + // Extract coupling params + auto &drag_pkg = pm->packages.Get("drag"); + const auto &fcp = drag_pkg->template Param("full_coupling_params"); + + // Packing and indexing + static auto desc = MakePackDescriptor( + resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + // Device copies of species properties + auto mu_s = fcp.mu_s; + auto sigma_s = fcp.sigma_s; + auto eps_s = fcp.eps_s; + auto dof_s = fcp.dof_s; + const Real kb_code = fcp.kb_code; + const GasDragModel cmodel = fcp.model; + + auto &artemis_pkg = pm->packages.Get("artemis"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "CoupleTwoFluids", DevExecSpace(), 0, md->NumBlocks() - 1, + kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + if (vmesh.GetSize(b, gas::cons::density()) < 2) return; + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + + // --- Extract species state --- + // Species 0 + Real &d0 = vmesh(b, gas::cons::density(0), k, j, i); + Real &e0 = vmesh(b, gas::cons::total_energy(0), k, j, i); + Real &u0 = vmesh(b, gas::cons::internal_energy(0), k, j, i); + d0 = std::max(d0, dflr_gas); + Real sie0 = ArtemisUtils::DualEnergySIE(vmesh, b, 0, k, j, i, de_switch, hx); + sie0 = std::max(sie0, sieflr_gas); + const Real T0_old = eos_d(0).TemperatureFromDensityInternalEnergy(d0, sie0); + const Real cv0_old = eos_d(0).SpecificHeatFromDensityTemperature(d0, T0_old); + + // Species 1 + Real &d1 = vmesh(b, gas::cons::density(1), k, j, i); + Real &e1 = vmesh(b, gas::cons::total_energy(1), k, j, i); + Real &u1 = vmesh(b, gas::cons::internal_energy(1), k, j, i); + d1 = std::max(d1, dflr_gas); + Real sie1 = ArtemisUtils::DualEnergySIE(vmesh, b, 1, k, j, i, de_switch, hx); + sie1 = std::max(sie1, sieflr_gas); + const Real T1_old = eos_d(1).TemperatureFromDensityInternalEnergy(d1, sie1); + const Real cv1_old = eos_d(1).SpecificHeatFromDensityTemperature(d1, T1_old); + + // --- Momentum coupling --- + // Compute Chapman-Cowing drag coefficient K_01 [mass/(vol*time)] + const Real T_pair = 0.5 * (T0_old + T1_old); + const Real K01 = + CollisionIntegrals::DragCoeff(cmodel, kb_code, mu_s(0), mu_s(1), sigma_s(0), + sigma_s(1), eps_s(0), eps_s(1), d0, d1, T_pair); + + const std::array v0{ + vmesh(b, gas::cons::momentum(VI(0, 0)), k, j, i) / (hx[0] * d0), + vmesh(b, gas::cons::momentum(VI(0, 1)), k, j, i) / (hx[1] * d0), + vmesh(b, gas::cons::momentum(VI(0, 2)), k, j, i) / (hx[2] * d0)}; + const std::array v1{ + vmesh(b, gas::cons::momentum(VI(1, 0)), k, j, i) / (hx[0] * d0), + vmesh(b, gas::cons::momentum(VI(1, 1)), k, j, i) / (hx[1] * d0), + vmesh(b, gas::cons::momentum(VI(1, 2)), k, j, i) / (hx[2] * d0)}; + + const Real ke0_old = 0.5 * d0 * (SQR(v0[0]) + SQR(v0[1]) + SQR(v0[2])); + const Real ke1_old = 0.5 * d1 * (SQR(v1[0]) + SQR(v1[1]) + SQR(v1[2])); + + for (int d = 0; d < 3; d++) { + const Real impulse = K01 * hx[d] * dt * (v1[d] - v0[d]) / + (1.0 + K01 * dt * (1.0 / d0 + 1.0 / d1)); + vmesh(b, gas::cons::momentum(VI(0, d)), k, j, i) += impulse; + vmesh(b, gas::cons::momentum(VI(1, d)), k, j, i) -= impulse; + } + const Real ke0_new = 0.5 * d0 * (SQR(v0[0]) + SQR(v0[1]) + SQR(v0[2])); + const Real ke1_new = 0.5 * d1 * (SQR(v1[0]) + SQR(v1[1]) + SQR(v1[2])); + + const Real dke0 = ke0_new - ke0_old; + const Real dke1 = ke1_new - ke1_old; + e0 += dke0; + e1 += dke1; + + const Real q_drag = + ClampRoundoffHeat(-(dke0 + dke1), std::abs(ke0_old) + std::abs(ke1_old)); + const Real c0_old = d0 * cv0_old; + const Real c1_old = d1 * cv1_old; + const Real csum_old = c0_old + c1_old + Fuzz(); + const Real q0_drag = q_drag * c0_old / csum_old; + const Real q1_drag = q_drag * c1_old / csum_old; + e0 += q0_drag; + e1 += q1_drag; + u0 += q0_drag; + u1 += q1_drag; + + // --- Thermal (energy) coupling --- + const Real sie0_post = std::max(u0 / d0, sieflr_gas); + const Real sie1_post = std::max(u1 / d1, sieflr_gas); + const Real T0 = eos_d(0).TemperatureFromDensityInternalEnergy(d0, sie0_post); + const Real T1 = eos_d(1).TemperatureFromDensityInternalEnergy(d1, sie1_post); + const Real cv0 = eos_d(0).SpecificHeatFromDensityTemperature(d0, T0); + const Real cv1 = eos_d(1).SpecificHeatFromDensityTemperature(d1, T1); + const Real H01 = CollisionIntegrals::ThermalRelaxRate( + cmodel, kb_code, mu_s(0), mu_s(1), sigma_s(0), sigma_s(1), eps_s(0), eps_s(1), + dof_s(0), dof_s(1), cv0, cv1, d0, d1, 0.5 * (T0 + T1)); + const Real c0 = d0 * cv0; + const Real c1 = d1 * cv1; + const Real q_th = + H01 * dt * (T1 - T0) / + (1.0 + H01 * dt * (1.0 / (c0 + Fuzz()) + 1.0 / (c1 + Fuzz()))); + e0 += q_th; + e1 -= q_th; + u0 += q_th; + u1 -= q_th; + + ApplyEnergyFloor(d0, sieflr_gas, ke0_new, e0, u0); + ApplyEnergyFloor(d1, sieflr_gas, ke1_new, e1, u1); + }); + return TaskStatus::complete; +} + +template +TaskStatus CoupleNFluids(MeshData *md, const int nmax, const Real dt) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + using TE = parthenon::TopologicalElement; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Extract gas package and params + auto &gas_pkg = pm->packages.Get("gas"); + const auto eos_d = gas_pkg->template Param>("eos_d"); + const auto dflr_gas = gas_pkg->template Param("dfloor"); + const auto sieflr_gas = gas_pkg->template Param("siefloor"); + const auto de_switch = gas_pkg->template Param("de_switch"); + const auto scr_level = gas_pkg->template Param("scr_level"); + + // Extract coupling params + auto &drag_pkg = pm->packages.Get("drag"); + const auto &fcp = drag_pkg->template Param("full_coupling_params"); + auto mu_s = fcp.mu_s; + auto sigma_s = fcp.sigma_s; + auto eps_s = fcp.eps_s; + auto dof_s = fcp.dof_s; + const Real kb_code = fcp.kb_code; + const GasDragModel cmodel = fcp.model; + + auto &artemis_pkg = pm->packages.Get("artemis"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + + // Packing and indexing + static auto desc = MakePackDescriptor( + resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + const int il = ib.s; + const int iu = ib.e; + + const int ncells1 = (iu - il + 1) + 2 * parthenon::Globals::nghost; + + // Scratch: velocity/temperature per species (nmax each), plus LU system (nmax x nmax) + // and a saved pre-coupling kinetic energy per species. + const int scr_size = + ScratchPad2D::shmem_size(ncells1, nmax) * 5 // vel, T, rho, rhs, ke_old + + ScratchPad3D::shmem_size(ncells1, nmax, nmax) // A + + ScratchPad2D::shmem_size(ncells1, nmax); // IPIV + + parthenon::par_for_outer( + DEFAULT_OUTER_LOOP_PATTERN, "CoupleNFluids", DevExecSpace(), scr_size, scr_level, 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, + KOKKOS_LAMBDA(parthenon::team_mbr_t mbr, const int b, const int k, const int j) { + const int ns = vmesh.GetSize(b, gas::cons::density()); + if (ns <= 1) return; + + // Allocate scratch + ScratchPad3D vel_s(mbr.team_scratch(scr_level), ncells1, 3, nmax); + ScratchPad2D T_s(mbr.team_scratch(scr_level), ncells1, nmax); + ScratchPad2D rho_s(mbr.team_scratch(scr_level), ncells1, nmax); + ScratchPad2D rhs(mbr.team_scratch(scr_level), ncells1, nmax); + ScratchPad2D ke_old_s(mbr.team_scratch(scr_level), ncells1, nmax); + ScratchPad3D A(mbr.team_scratch(scr_level), ncells1, nmax, nmax); + ScratchPad2D IPIV(mbr.team_scratch(scr_level), ncells1, nmax); + + // ----------------------------------------------------------------- + // Fill: load species state and assemble momentum coupling system. + // We solve each spatial direction independently with the same K_ij matrix. + // The system is: + // (I + dt * K_hat) * v_new = v_old + // where K_hat_{nn} = sum_{m!=n} K_nm / rho_n + // K_hat_{nm} = -K_nm / rho_n (m != n) + // This is solved direction-by-direction by reusing A (same matrix per direction). + // ----------------------------------------------------------------- + + // Step 1: load densities, temperatures, floor + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + for (int n = 0; n < ns; ++n) { + Real &dens = vmesh(b, gas::cons::density(n), k, j, i); + dens = std::max(dens, dflr_gas); + rho_s(i, n) = dens; + Real sie = + ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + sie = std::max(sie, sieflr_gas); + T_s(i, n) = eos_d(n).TemperatureFromDensityInternalEnergy(dens, sie); + Real ke = 0.0; + for (int d = 0; d < 3; d++) { + const Real v = + vmesh(b, gas::cons::momentum(VI(n, d)), k, j, i) / (hx[d] * dens); + vel_s(i, d, n) = v; + ke += SQR(v); + } + ke_old_s(i, n) = 0.5 * dens * ke; + } + }); + mbr.team_barrier(); + + // Step 2: build coupling matrix A and solve for each momentum direction + for (int dir = 0; dir < 3; ++dir) { + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + // Build A (identity + dt * K_hat) and rhs (= v_old) + for (int n = 0; n < ns; ++n) { + const Real &rho_n = rho_s(i, n); + const Real &T_n = T_s(i, n); + // Diagonal: 1 + sum_m K_nm / rho_n + Real diag = 1.0; + for (int m = 0; m < ns; ++m) { + if (m == n) continue; + const Real T_pair = 0.5 * (T_n + T_s(i, m)); + const Real K_nm = CollisionIntegrals::DragCoeff( + cmodel, kb_code, mu_s(n), mu_s(m), sigma_s(n), sigma_s(m), + eps_s(n), eps_s(m), rho_n, rho_s(i, m), T_pair); + diag += dt * K_nm / rho_n; + A(i, n, m) = -dt * K_nm / rho_n; + } + A(i, n, n) = diag; + IPIV(i, n) = 0; + rhs(i, n) = vel_s(i, dir, n); + } + }); + mbr.team_barrier(); + + // Solve A * v_new = v_old + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + auto A_ = Kokkos::subview(A, i, Kokkos::ALL, Kokkos::ALL); + auto IPIV_ = Kokkos::subview(IPIV, i, Kokkos::ALL); + auto RHS_ = Kokkos::subview(rhs, i, Kokkos::ALL); + KokkosBatched::SerialGetrf::invoke( + A_, IPIV_); + KokkosBatched::SerialGetrs< + KokkosBatched::Trans::NoTranspose, + KokkosBatched::Algo::Getrs::Unblocked>::invoke(A_, IPIV_, RHS_); + }); + mbr.team_barrier(); + + // Apply momentum updates + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + for (int n = 0; n < ns; ++n) { + const Real rho_n = rho_s(i, n); + const Real v_old = vel_s(i, dir, n); + const Real v_new = rhs(i, n); // solution overwrites rhs + const Real dv = v_new - v_old; + vmesh(b, gas::cons::momentum(VI(n, dir)), k, j, i) += + rho_n * hx[dir] * dv; + } + }); + mbr.team_barrier(); + } // dir loop + + // Deposit kinetic energy changes and irreversible drift heating. + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + Real q_drag = 0.0; + Real ke_scale = 0.0; + Real csum = 0.0; + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + for (int n = 0; n < ns; ++n) { + const Real &rho_n = rho_s(i, n); + const Real cv_n = + eos_d(n).SpecificHeatFromDensityTemperature(rho_n, T_s(i, n)); + const std::array v0{ + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * rho_n), + vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * rho_n), + vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * rho_n)}; + + const Real ke_new = 0.5 * rho_n * (SQR(v0[0]) + SQR(v0[1]) + SQR(v0[2])); + const Real dke = ke_new - ke_old_s(i, n); + vmesh(b, gas::cons::total_energy(n), k, j, i) += dke; + q_drag -= dke; + ke_scale += std::abs(ke_old_s(i, n)); + rhs(i, n) = rho_n * cv_n; + csum += rhs(i, n); + } + q_drag = ClampRoundoffHeat(q_drag, ke_scale); + for (int n = 0; n < ns; ++n) { + const Real qn = q_drag * rhs(i, n) / (csum + Fuzz()); + vmesh(b, gas::cons::total_energy(n), k, j, i) += qn; + vmesh(b, gas::cons::internal_energy(n), k, j, i) += qn; + } + }); + mbr.team_barrier(); + + // Refresh temperatures after drift heating before thermal relaxation. + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + for (int n = 0; n < ns; ++n) { + const Real rho_n = rho_s(i, n); + const Real sie = std::max( + vmesh(b, gas::cons::internal_energy(n), k, j, i) / rho_n, sieflr_gas); + T_s(i, n) = eos_d(n).TemperatureFromDensityInternalEnergy(rho_n, sie); + } + }); + mbr.team_barrier(); + + // Thermal relaxation: same impicit pattern with scalar exchange matrix. + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + for (int n = 0; n < ns; ++n) { + const Real rho_n = rho_s(i, n); + const Real T_n = T_s(i, n); + const Real cv_n = eos_d(n).SpecificHeatFromDensityTemperature(rho_n, T_n); + rhs(i, n) = T_n; + Real diag = 1.0; + for (int m = 0; m < ns; ++m) { + if (m == n) continue; + const Real rho_m = rho_s(i, m); + const Real T_m = T_s(i, m); + const Real cv_m = + eos_d(m).SpecificHeatFromDensityTemperature(rho_m, T_m); + const Real H_nm = CollisionIntegrals::ThermalRelaxRate( + cmodel, kb_code, mu_s(n), mu_s(m), sigma_s(n), sigma_s(m), eps_s(n), + eps_s(m), dof_s(n), dof_s(m), cv_n, cv_m, rho_n, rho_m, + 0.5 * (T_n + T_m)); + const Real rate = dt * H_nm / (rho_n * cv_n + Fuzz()); + diag += rate; + A(i, n, m) = -rate; + } + A(i, n, n) = diag; + IPIV(i, n) = 0; + } + }); + mbr.team_barrier(); + + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + auto A_ = Kokkos::subview(A, i, Kokkos::ALL, Kokkos::ALL); + auto IPIV_ = Kokkos::subview(IPIV, i, Kokkos::ALL); + auto RHS_ = Kokkos::subview(rhs, i, Kokkos::ALL); + KokkosBatched::SerialGetrf::invoke( + A_, IPIV_); + KokkosBatched::SerialGetrs< + KokkosBatched::Trans::NoTranspose, + KokkosBatched::Algo::Getrs::Unblocked>::invoke(A_, IPIV_, RHS_); + }); + mbr.team_barrier(); + + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, mbr, il, iu, [&](const int i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + for (int n = 0; n < ns; ++n) { + const Real rho_n = rho_s(i, n); + const Real T_old = T_s(i, n); + const Real T_new = rhs(i, n); + const Real cv_n = + eos_d(n).SpecificHeatFromDensityTemperature(rho_n, T_old); + const Real dE = rho_n * cv_n * (T_new - T_old); + vmesh(b, gas::cons::total_energy(n), k, j, i) += dE; + vmesh(b, gas::cons::internal_energy(n), k, j, i) += dE; + const std::array v0{ + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) / (hx[0] * rho_n), + vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) / (hx[1] * rho_n), + vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) / (hx[2] * rho_n)}; + const Real ke_new = 0.5 * rho_n * (SQR(v0[0]) + SQR(v0[1]) + SQR(v0[2])); + + Real &etot = vmesh(b, gas::cons::total_energy(n), k, j, i); + Real &eint = vmesh(b, gas::cons::internal_energy(n), k, j, i); + ApplyEnergyFloor(rho_n, sieflr_gas, ke_new, etot, eint); + } + }); + }); // par_for_outer + return TaskStatus::complete; +} + +} // namespace Drag + +#endif // DRAG_DRAG_IMPL_HPP_ diff --git a/src/dust/coagulation/coagulation.cpp b/src/dust/coagulation/coagulation.cpp index 95dd7f56..532b63b3 100644 --- a/src/dust/coagulation/coagulation.cpp +++ b/src/dust/coagulation/coagulation.cpp @@ -222,7 +222,7 @@ TaskStatus CoagulationStep(MeshData *md, const Real time, const Real dt) { // Extract EOS auto &gas_pkg = pm->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Extract dust params auto &dust_pkg = pm->packages.Get("dust"); @@ -306,8 +306,8 @@ TaskStatus CoagulationStep(MeshData *md, const Real time, const Real dt) { // Extract gas state vector const Real &gdens = vmesh(b, gas::prim::density(0), k, j, i); const Real &gsie = vmesh(b, gas::prim::sie(0), k, j, i); - const Real kT = eos_d.TemperatureFromDensityInternalEnergy(gdens, gsie); - const Real &gbulk = eos_d.BulkModulusFromDensityInternalEnergy(gdens, gsie); + const Real kT = eos_d(0).TemperatureFromDensityInternalEnergy(gdens, gsie); + const Real &gbulk = eos_d(0).BulkModulusFromDensityInternalEnergy(gdens, gsie); const Real cs1 = std::sqrt(gbulk / gdens) * vel0; const Real gdens1 = gdens * rho0; const Real kT1 = kT * kT0; diff --git a/src/gas/cooling/beta_cooling.cpp b/src/gas/cooling/beta_cooling.cpp index b1aa6e36..fa588c3c 100644 --- a/src/gas/cooling/beta_cooling.cpp +++ b/src/gas/cooling/beta_cooling.cpp @@ -47,7 +47,7 @@ TaskStatus BetaCooling(MeshData *md, const Real time, const Real dt) { // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); - const auto &eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); const auto de_switch = gas_pkg->template Param("de_switch"); const auto dflr_gas = gas_pkg->template Param("dfloor"); const auto sieflr_gas = gas_pkg->template Param("siefloor"); @@ -131,8 +131,8 @@ TaskStatus BetaCooling(MeshData *md, const Real time, const Real dt) { sie = (efloor)*sie + (!efloor) * sieflr_gas; // Compute the energy change from the temperature change - const Real cv = eos_d.SpecificHeatFromDensityInternalEnergy(dens, sie); - const Real Tn = eos_d.TemperatureFromDensityInternalEnergy(dens, sie); + const Real cv = eos_d(n).SpecificHeatFromDensityInternalEnergy(dens, sie); + const Real Tn = eos_d(n).TemperatureFromDensityInternalEnergy(dens, sie); const Real dE = -dens * cv * omdt / (beta + omdt) * (Tn - T0); // Add this energy change to the conserved fields diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index f643b472..aa84c3b0 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -48,6 +48,13 @@ std::shared_ptr Initialize(ParameterInput *pin, auto gas = std::make_shared("gas"); Params ¶ms = gas->AllParams(); + // Number of gas species + const int nspecies = pin->GetOrAddInteger("gas", "nspecies", 1); + params.Add("nspecies", nspecies); + std::vector fluidids; + for (int n = 0; n < nspecies; ++n) + fluidids.push_back(n); + // Fluid behavior for this package Fluid fluid_type = Fluid::gas; params.Add("fluid_type", fluid_type); @@ -79,35 +86,71 @@ std::shared_ptr Initialize(ParameterInput *pin, if (pin->DoesBlockExist("gas/eos/ideal") || (pin->DoesParameterExist("gas", "gamma"))) { const std::string block_name = pin->DoesBlockExist("gas/eos/ideal") ? "gas/eos/ideal" : "gas"; - const Real gamma = pin->GetOrAddReal(block_name, "gamma", 1.66666666667); - auto cv = Null(); - auto mu = Null(); + std::vector gamma_default(nspecies, 1.66666666667); + std::vector mu_default(nspecies, Null()); + std::vector cv_default(nspecies, Null()); + + auto gamma_v = pin->GetOrAddVector(block_name, "gamma", gamma_default); + PARTHENON_REQUIRE(gamma_v.size() == static_cast(nspecies), + "gamma must have nspecies entries"); + auto cv_v = cv_default; + auto mu_v = mu_default; if (pin->DoesParameterExist(block_name, "cv")) { PARTHENON_REQUIRE(!pin->DoesParameterExist("gas", "mu"), "Cannot specify both cv and mu"); - cv = pin->GetReal(block_name, "cv"); - PARTHENON_REQUIRE(cv > 0, "Only positive cv allowed!"); - mu = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * cv); + cv_v = pin->GetVector(block_name, "cv"); + PARTHENON_REQUIRE(cv_v.size() == static_cast(nspecies), + "cv must have nspecies entries"); + for (int n = 0; n < nspecies; ++n) { + PARTHENON_REQUIRE(cv_v[n] > 0, "Only positive cv allowed!"); + mu_v[n] = constants.GetKBCode() / + ((gamma_v[n] - 1.) * constants.GetAMUCode() * cv_v[n]); + } } else { - mu = pin->GetOrAddReal(block_name, "mu", 1.); - PARTHENON_REQUIRE(mu > 0, "Only positive mean molecular weight allowed!"); - cv = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * mu); + std::vector mu_ones(nspecies, 1.); + mu_v = pin->GetOrAddVector(block_name, "mu", mu_ones); + PARTHENON_REQUIRE(mu_v.size() == static_cast(nspecies), + "mu must have nspecies entries"); + for (int n = 0; n < nspecies; ++n) { + PARTHENON_REQUIRE(mu_v[n] > 0, "Only positive mean molecular weight allowed!"); + cv_v[n] = constants.GetKBCode() / + ((gamma_v[n] - 1.) * constants.GetAMUCode() * mu_v[n]); + } + } + + ParArray1D gamma("gamma", nspecies); + ParArray1D mu("mu", nspecies); + ParArray1D cv("cv", nspecies); + auto h_gamma = gamma.GetHostMirror(); + auto h_mu = mu.GetHostMirror(); + auto h_cv = cv.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + h_gamma(n) = gamma_v[n]; + h_mu(n) = mu_v[n]; + h_cv(n) = cv_v[n]; } + gamma.DeepCopy(h_gamma); + mu.DeepCopy(h_mu); + cv.DeepCopy(h_cv); + eos_type = "ideal"; - params.Add("kbmu", constants.GetKBCode() / (mu * constants.GetAMUCode())); + params.Add("kbmu", constants.GetKBCode() / (mu_v[0] * constants.GetAMUCode())); params.Add("mu", mu); params.Add("cv", cv); params.Add("kb", constants.GetKBCode()); params.Add("amu", constants.GetAMUCode()); - params.Add("Rgas", constants.GetKBCode() / (constants.GetAMUCode() * mu)); - - // Build EOS - EOS eos_host = singularity::UnitSystem( - singularity::IdealGas(gamma - 1., cv * units.GetSpecificHeatCodeToPhysical()), - singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimePhysicalToCode(), - units.GetMassPhysicalToCode(), units.GetLengthPhysicalToCode(), - units.GetTemperaturePhysicalToCode()); - EOS eos_device = eos_host.GetOnDevice(); + params.Add("Rgas", constants.GetKBCode() / (constants.GetAMUCode() * mu_v[0])); + ParArray1D eos_device("eos_d", nspecies); + auto eos_host = eos_device.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + eos_host(n) = singularity::UnitSystem( + singularity::IdealGas(gamma_v[n] - 1., + cv_v[n] * units.GetSpecificHeatCodeToPhysical()), + singularity::eos_units_init::LengthTimeUnitsInit(), + units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), + units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); + eos_device(n) = eos_host(n).GetOnDevice(); + } params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); // TODO This needs to be removed when we convert everything to EOS calls @@ -118,76 +161,145 @@ std::shared_ptr Initialize(ParameterInput *pin, const std::string block_name = "gas/eos/h-he"; if (pin->DoesParameterExist(block_name, "eos_file")) { // load from file - const std::string filename = pin->GetString(block_name, "eos_file"); - EOS eos_host = singularity::UnitSystem( - ArtemisEOS::IdealHHe(filename), - singularity::eos_units_init::LengthTimeUnitsInit(), - units.GetTimePhysicalToCode(), units.GetMassPhysicalToCode(), - units.GetLengthPhysicalToCode(), units.GetTemperaturePhysicalToCode()); - EOS eos_device = eos_host.GetOnDevice(); - params.Add("mu", pin->GetOrAddReal(block_name, "mu", 1.)); + std::vector filenames(nspecies, ""); + pin->GetOrAddVector(block_name, "eos_file", filenames); + ParArray1D eos_device("eos_d", nspecies); + ParArray1D eos_host = eos_device.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + eos_host(n) = singularity::UnitSystem( + ArtemisEOS::IdealHHe(filenames[n]), + singularity::eos_units_init::LengthTimeUnitsInit(), + units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), + units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); + eos_device(n) = eos_host(n).GetOnDevice(); + } + std::vector mu_default(nspecies, 1.); + auto mu_v = pin->GetOrAddVector(block_name, "mu", mu_default); + PARTHENON_REQUIRE(mu_v.size() == static_cast(nspecies), + "mu must have nspecies entries"); + ParArray1D mu("mu", nspecies); + auto h_mu = mu.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + h_mu(n) = mu_v[n]; + } + mu.DeepCopy(h_mu); + params.Add("mu", mu); params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); } else { const std::string save_to_file = pin->GetOrAddString(block_name, "save_to_file", ""); - const Real X = pin->GetReal(block_name, "x"); - const Real Y = pin->GetReal(block_name, "y"); const Real ltmin = pin->GetOrAddReal(block_name, "ltmin", 0); const Real ltmax = pin->GetOrAddReal(block_name, "ltmax", 6); const Real ldmin = pin->GetOrAddReal(block_name, "ldmin", -15); const Real ldmax = pin->GetOrAddReal(block_name, "ldmax", -3); const int nd = pin->GetOrAddInteger(block_name, "nd", 100); const int nt = pin->GetOrAddInteger(block_name, "nt", 100); - ArtemisEOS::IdealHHe eos_base(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, - save_to_file, true); - eos_base.SetFloors(siefloor, 0.0, dfloor, 0.0); - - EOS eos_host = singularity::UnitSystem( - std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), - units.GetTimePhysicalToCode(), units.GetMassPhysicalToCode(), - units.GetLengthPhysicalToCode(), units.GetTemperaturePhysicalToCode()); - EOS eos_device = eos_host.GetOnDevice(); - params.Add("mu", pin->GetOrAddReal(block_name, "mu", 1.)); + + std::vector X_v(nspecies, 0.); + std::vector Y_v(nspecies, 0.); + X_v = pin->GetVector(block_name, "x"); + Y_v = pin->GetVector(block_name, "y"); + PARTHENON_REQUIRE(X_v.size() == static_cast(nspecies), + "x must have nspecies entries"); + PARTHENON_REQUIRE(Y_v.size() == static_cast(nspecies), + "y must have nspecies entries"); + + ParArray1D eos_device("eos_d", nspecies); + auto eos_host = eos_device.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + ArtemisEOS::IdealHHe eos_base(X_v[n], Y_v[n], ltmin, ltmax, nt, ldmin, ldmax, nd, + save_to_file, true); + eos_base.SetFloors(siefloor, 0.0, dfloor, 0.0); + + eos_host(n) = singularity::UnitSystem( + std::move(eos_base), singularity::eos_units_init::LengthTimeUnitsInit(), + units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), + units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); + eos_device(n) = eos_host(n).GetOnDevice(); + } + auto mu_v = + pin->GetOrAddVector(block_name, "mu", std::vector(nspecies, 1.)); + PARTHENON_REQUIRE(mu_v.size() == static_cast(nspecies), + "mu must have nspecies entries"); + ParArray1D mu("mu", nspecies); + auto h_mu = mu.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + h_mu(n) = mu_v[n]; + } + mu.DeepCopy(h_mu); + params.Add("mu", mu); params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); } -#ifdef WITH_SESAME - } else if (pin->DoesBlockExist("gas/eos/sesame_re")) { - eos_type = "sesame_re"; + } else if (pin->DoesBlockExist("gas/eos/spiner_re")) { + eos_type = "spiner_re"; params.Add("eos_type", eos_type); - const std::string block_name = "gas/eos/sesame_re"; - std::string filename = pin->GetString(block_name, "eos_file"); - EOS eos_host = singularity::UnitSystem( - singularity::SpinerEOSDependsRhoSie(filename, "gas"), - singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimePhysicalToCode(), - units.GetMassPhysicalToCode(), units.GetLengthPhysicalToCode(), - units.GetTemperaturePhysicalToCode()); - EOS eos_device = eos_host.GetOnDevice(); + const std::string block_name = "gas/eos/spiner_re"; + auto filenames = pin->GetVector(block_name, "eos_file"); + auto matid = pin->GetVector(block_name, "matid"); + PARTHENON_REQUIRE(filenames.size() == static_cast(nspecies), + "eos_file must have nspecies entries"); + ParArray1D eos_device("eos_d", nspecies); + auto eos_host = eos_device.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + eos_host(n) = singularity::UnitSystem( + singularity::SpinerEOSDependsRhoSie(filenames[n], matid[n]), + singularity::eos_units_init::LengthTimeUnitsInit(), + units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), + units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); + eos_device(n) = eos_host(n).GetOnDevice(); + } + std::vector mu_default(nspecies, 1.); + auto mu_v = pin->GetOrAddVector(block_name, "mu", mu_default); + PARTHENON_REQUIRE(mu_v.size() == static_cast(nspecies), + "mu must have nspecies entries"); + ParArray1D mu("mu", nspecies); + auto h_mu = mu.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + h_mu(n) = mu_v[n]; + } + mu.DeepCopy(h_mu); + params.Add("mu", mu); params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); - params.Add("mu", pin->GetOrAddReal(block_name, "mu", 1.)); - } else if (pin->DoesBlockExist("gas/eos/sesame_rt")) { - eos_type = "sesame_rt"; - const std::string block_name = "gas/eos/sesame_rt"; - std::string filename = pin->GetString(block_name, "eos_file"); - EOS eos_host = singularity::UnitSystem( - singularity::SpinerEOSDependsRhoT(filename, "gas"), - singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimePhysicalToCode(), - units.GetMassPhysicalToCode(), units.GetLengthPhysicalToCode(), - units.GetTemperaturePhysicalToCode()); - EOS eos_device = eos_host.GetOnDevice(); + } else if (pin->DoesBlockExist("gas/eos/spiner_rt")) { + eos_type = "spiner_rt"; + const std::string block_name = "gas/eos/spiner_rt"; + auto filenames = pin->GetVector(block_name, "eos_file"); + auto matid = pin->GetVector(block_name, "matid"); + + PARTHENON_REQUIRE(filenames.size() == static_cast(nspecies), + "eos_file must have nspecies entries"); + ParArray1D eos_device("eos_d", nspecies); + auto eos_host = eos_device.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + eos_host(n) = singularity::UnitSystem( + singularity::SpinerEOSDependsRhoT(filenames[n], matid[n]), + singularity::eos_units_init::LengthTimeUnitsInit(), + units.GetTimeCodeToPhysical(), units.GetMassCodeToPhysical(), + units.GetLengthCodeToPhysical(), units.GetTemperatureCodeToPhysical()); + eos_device(n) = eos_host(n).GetOnDevice(); + } + std::vector mu_default(nspecies, 1.); + auto mu_v = pin->GetOrAddVector(block_name, "mu", mu_default); + PARTHENON_REQUIRE(mu_v.size() == static_cast(nspecies), + "mu must have nspecies entries"); + ParArray1D mu("mu", nspecies); + auto h_mu = mu.GetHostMirror(); + for (int n = 0; n < nspecies; ++n) { + h_mu(n) = mu_v[n]; + } + mu.DeepCopy(h_mu); + params.Add("mu", mu); params.Add("eos_h", eos_host); params.Add("eos_d", eos_device); - params.Add("mu", pin->GetOrAddReal(block_name, "mu", 1.)); -#endif // WITH_SESAME #endif // SPINER_USE_HDF } else { PARTHENON_FAIL("Unsupported gas EOS!"); } params.Add("eos_type", eos_type); - // Riemann solver RSolver riemann_solver = RSolver::null; const std::string riemann = pin->GetOrAddString("gas", "riemann", "hllc-general"); @@ -326,21 +438,29 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("cond_params", dp); } - // Number of gas species - const int nspecies = pin->GetOrAddInteger("gas", "nspecies", 1); - params.Add("nspecies", nspecies); - std::vector fluidids; - for (int n = 0; n < nspecies; ++n) - fluidids.push_back(n); - // Scratch for gas flux const int scr_level = pin->GetOrAddInteger("gas", "scr_level", 0); params.Add("scr_level", scr_level); + params.Add("radius", pin->GetOrAddInteger("gas", "radius_cgs", 0.0) * + units.GetLengthPhysicalToCode()); + // Logarithmic gridding? const bool log = pin->GetOrAddString("artemis", "radial_spacing", "uniform") == "logarithmic"; + Real sparse_alloc_thresh = + pin->GetOrAddReal("gas", "density_allocation_threshold", 1e-30); + Real sparse_dealloc_thresh = + pin->GetOrAddReal("gas", "density_deallocation_threshold", 1e-32); + if (nspecies == 1) { + sparse_alloc_thresh = 0.0; + sparse_dealloc_thresh = 0.0; + } + + params.Add("sparse_alloc_thresh", sparse_alloc_thresh); + params.Add("sparse_dealloc_thresh", sparse_dealloc_thresh); + // Control field for sparse gas fields std::string control_field = gas::cons::density::name(); @@ -348,7 +468,7 @@ std::shared_ptr Initialize(ParameterInput *pin, Metadata m = Metadata({Metadata::Cell, Metadata::Conserved, Metadata::Independent, Metadata::WithFluxes, Metadata::Sparse}); ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); - m.SetSparseThresholds(0.0, 0.0, 0.0); + m.SetSparseThresholds(sparse_alloc_thresh, sparse_dealloc_thresh, 0.0); gas->AddSparsePool(m, control_field, fluidids); // Conserved Momenta @@ -542,7 +662,7 @@ Real EstimateTimestepMesh(MeshData *md) { auto &gas_pkg = pm->packages.Get("gas"); auto ¶ms = gas_pkg->AllParams(); - auto eos_d = params.template Get("eos_d"); + const auto &eos_d = params.template Get>("eos_d"); static auto desc = MakePackDescriptor *md); TaskStatus CalculateFluxes(MeshData *md, const bool pcm); TaskStatus FluxSource(MeshData *md, const Real dt); +TaskStatus ApplyClosure(MeshData *md, const Real dt); template TaskStatus DiffusionUpdate(MeshData *md, const Real dt); diff --git a/src/pgen/beam.hpp b/src/pgen/beam.hpp index 5a1d71fb..27b1bc0e 100644 --- a/src/pgen/beam.hpp +++ b/src/pgen/beam.hpp @@ -67,7 +67,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool do_rad = artemis_pkg->Param("do_moment"); auto rad_pkg = pmb->packages.Get("moments"); auto gas_pkg = pmb->packages.Get("gas"); - const auto eos = gas_pkg->Param("eos_d"); + const auto eos = gas_pkg->Param>("eos_d"); const Real ar = rad_pkg->Param("arad"); // packing and capture variables for kernel @@ -106,7 +106,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(0, gas::prim::velocity(1), k, j, i) = 0.0; v(0, gas::prim::velocity(2), k, j, i) = 0.0; v(0, gas::prim::sie(0), k, j, i) = - eos.InternalEnergyFromDensityTemperature(dens, pars.tg); + eos(0).InternalEnergyFromDensityTemperature(dens, pars.tg); } if (do_rad) { v(0, rad::prim::energy(0), k, j, i) = ar * SQR(SQR(pars.tg)); diff --git a/src/pgen/blast.hpp b/src/pgen/blast.hpp index bba557c8..3d8c0853 100644 --- a/src/pgen/blast.hpp +++ b/src/pgen/blast.hpp @@ -154,7 +154,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool do_dust = artemis_pkg->Param("do_dust"); // TODO(PDM): Replace the below with a call to singularity-eos auto gas_pkg = pmb->packages.Get("gas"); - const auto &eos = gas_pkg->Param("eos_d"); + const auto &eos = gas_pkg->template Param>("eos_d"); // packing and capture variables for kernel auto &md = pmb->meshblock_data.Get(); @@ -184,7 +184,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real total_vol = coords.GetVolume(vg, 0, k, j, i); const auto &xv = coords.GetCellCenter(vg, 0, k, j, i); Real den = pars.d0; - Real e0 = ArtemisUtils::EofPR(eos, pars.p0, den); + Real e0 = ArtemisUtils::EofPR(eos(0), pars.p0, den); Real internal_energy = 0.0; auto xcart = coords.ConvertToCart(xv); const auto &xc = coords.ConvertToCart(pars.x0); diff --git a/src/pgen/coag.hpp b/src/pgen/coag.hpp index a73a4a76..9507b7ac 100644 --- a/src/pgen/coag.hpp +++ b/src/pgen/coag.hpp @@ -91,7 +91,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { dcv.rho0 = pin->GetOrAddReal("problem", "rho0", 1.0); auto gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Extract adiabatic index and H0 dcv.gamma = gas_pkg->Param("adiabatic_index"); @@ -104,8 +104,8 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const auto &constants = artemis_pkg->Param("constants"); const Real kbmu = constants.GetKBCode() / (mu * constants.GetAMUCode()); const Real gtemp = SQR(dcv.h0) / kbmu / dcv.gamma; - const Real gsie = eos_d.InternalEnergyFromDensityTemperature(gdens, gtemp); - const Real pres = eos_d.PressureFromDensityTemperature(gdens, gtemp); + const Real gsie = eos_d(0).InternalEnergyFromDensityTemperature(gdens, gtemp); + const Real pres = eos_d(0).PressureFromDensityTemperature(gdens, gtemp); if (pmb->gid == 0) { std::cout << "gamma, h0, pres=" << dcv.gamma << " " << dcv.h0 << " " << dcv.gm1 * gdens * gsie << " " << pres << std::endl; diff --git a/src/pgen/conduction.hpp b/src/pgen/conduction.hpp index 91dc9189..335bcafb 100644 --- a/src/pgen/conduction.hpp +++ b/src/pgen/conduction.hpp @@ -69,13 +69,13 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool do_gas = artemis_pkg->Param("do_gas"); const bool do_gravity = artemis_pkg->Param("do_gravity"); auto cond_params = artemis_pkg->Param("cond_pgen_params"); - EOS eos_d; + ParArray1D eos_d; Real gx1 = 0.; if (do_gas) { auto gas_pkg = pmb->packages.Get("gas"); PARTHENON_REQUIRE(gas_pkg->Param("nspecies") == 1, "Cond pgen requires a single gas species.") - eos_d = gas_pkg->Param("eos_d"); + eos_d = gas_pkg->Param>("eos_d"); if (do_gravity) { auto grav_pkg = pmb->packages.Get("gravity"); if (grav_pkg->Param("type") == @@ -111,7 +111,8 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { geometry::Coords coords(cpars, pco, k, j, i); const auto &xv = coords.GetCellCenter(); - const Real P0 = eos_d.PressureFromDensityTemperature(pars.g_rho, pars.g_temp); + const Real P0 = + eos_d(0).PressureFromDensityTemperature(pars.g_rho, pars.g_temp); const Real Rgas = P0 / (pars.g_rho * pars.g_temp); const Real P = P0 * std::exp(gx1 * pars.g_rho / P0 * (xv[0] - x1min)); const Real dens = P / (Rgas * pars.g_temp); @@ -121,7 +122,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(0, gas::prim::velocity(1), k, j, i) = pars.g_vx2; v(0, gas::prim::velocity(2), k, j, i) = pars.g_vx3; v(0, gas::prim::sie(0), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(dens, pars.g_temp); + eos_d(0).InternalEnergyFromDensityTemperature(dens, pars.g_temp); } }); } @@ -142,7 +143,7 @@ void CondBoundaryImpl(std::shared_ptr> &mbd, bool coarse) { auto &gas_pkg = pmb->packages.Get("gas"); auto do_gas = artemis_pkg->template Param("do_gas"); auto diff_params = gas_pkg->Param("cond_params"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Packing static auto descriptors = @@ -231,10 +232,10 @@ void CondBoundaryImpl(std::shared_ptr> &mbd, bool coarse) { // active zone const Real da = v(0, gas::prim::density(0), ia[0], ia[1], ia[2]); const Real siea = v(0, gas::prim::sie(0), ia[0], ia[1], ia[2]); - const Real Ta = eos_d.TemperatureFromDensityInternalEnergy(da, siea); - const Real Pa = eos_d.PressureFromDensityInternalEnergy(da, siea); + const Real Ta = eos_d(0).TemperatureFromDensityInternalEnergy(da, siea); + const Real Pa = eos_d(0).PressureFromDensityInternalEnergy(da, siea); - const Real ka = dcoeff.Get(dcp, ca, xva, da, siea, eos_d); + const Real ka = dcoeff.Get(dcp, ca, xva, da, siea, eos_d(0)); Real Tg = dp.g_temp; if (INNER) { Tg = Ta - dp.flux * xma / ka; @@ -242,7 +243,7 @@ void CondBoundaryImpl(std::shared_ptr> &mbd, bool coarse) { // Density from dP/dx = - rho g const Real densg = da * (Ta - 0.5 * gx1 * xma) / (Tg + 0.5 * gx1 * xma); - const Real sieg = eos_d.InternalEnergyFromDensityTemperature(densg, Tg); + const Real sieg = eos_d(0).InternalEnergyFromDensityTemperature(densg, Tg); // Extrapolate gas velocity Real gva[3] = {v(0, gas::prim::velocity(0), ia[0], ia[1], ia[2]), diff --git a/src/pgen/constant.hpp b/src/pgen/constant.hpp index a7181a1a..52b950d8 100644 --- a/src/pgen/constant.hpp +++ b/src/pgen/constant.hpp @@ -59,12 +59,12 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_gas = artemis_pkg->Param("do_gas"); const bool do_dust = artemis_pkg->Param("do_dust"); - EOS eos_d; + ParArray1D eos_d; if (do_gas) { auto gas_pkg = pmb->packages.Get("gas"); PARTHENON_REQUIRE(gas_pkg->Param("nspecies") == 1, "Constant pgen requires a single gas species.") - eos_d = gas_pkg->Param("eos_d"); + eos_d = gas_pkg->Param>("eos_d"); constant_params.g_rho = pin->GetOrAddReal("problem", "gas_rho", 1.0); constant_params.g_vx1 = pin->GetOrAddReal("problem", "gas_vx1", 0.0); constant_params.g_vx2 = pin->GetOrAddReal("problem", "gas_vx2", 0.0); @@ -151,7 +151,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(0, gas::prim::velocity(2), k, j, i) = (pars.g_vx1 * ex3[0] + pars.g_vx2 * ex3[1] + pars.g_vx3 * ex3[2]); v(0, gas::prim::sie(0), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(pars.g_rho, pars.g_temp); + eos_d(0).InternalEnergyFromDensityTemperature(pars.g_rho, pars.g_temp); } if (do_dust) { for (int n = 0; n < v.GetSize(0, dust::prim::density()); ++n) { diff --git a/src/pgen/crooked_pipe.hpp b/src/pgen/crooked_pipe.hpp index 45868151..2e6d6826 100644 --- a/src/pgen/crooked_pipe.hpp +++ b/src/pgen/crooked_pipe.hpp @@ -50,7 +50,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(!(update_fluxes), "Crooked pipe problem requires update fluxes to be off!"); auto gas_pkg = pmb->packages.Get("gas"); - const auto eos = gas_pkg->Param("eos_d"); + const auto eos = gas_pkg->Param>("eos_d"); Real ar = Null(); if (do_moment) { auto rad_pkg = pmb->packages.Get("moments"); @@ -118,21 +118,21 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // default thin cell v(0, gas::prim::density(), k, j, i) = rho_thin; v(0, gas::prim::sie(), k, j, i) = - eos.InternalEnergyFromDensityTemperature(rho_thin, t_init); + eos(0).InternalEnergyFromDensityTemperature(rho_thin, t_init); for (const auto &iregion : thick_regions) { if (xl >= iregion[0] && xu <= iregion[1] && yl >= iregion[2] && yu <= iregion[3]) { v(0, gas::prim::density(), k, j, i) = rho_thick; v(0, gas::prim::sie(), k, j, i) = - eos.InternalEnergyFromDensityTemperature(rho_thick, t_init); + eos(0).InternalEnergyFromDensityTemperature(rho_thick, t_init); } } for (const auto &iregion : thin_source_regions) { if (xl >= iregion[0] && xu <= iregion[1] && yl >= iregion[2] && yu <= iregion[3]) { v(0, gas::prim::sie(), k, j, i) = - eos.InternalEnergyFromDensityTemperature(rho_thin, t_source); + eos(0).InternalEnergyFromDensityTemperature(rho_thin, t_source); } } }); diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index fd089728..4d0ce399 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -277,7 +277,8 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { disk_params.temp_soft2 = pin->GetOrAddReal("problem", "temp_soft", 0.0); const auto mu = gas_pkg->Param("mu"); auto &constants = artemis_pkg->Param("constants"); - const auto &eos = gas_pkg->Param("eos_h"); + const auto &eos_arr = gas_pkg->Param>("eos_h"); + const auto &eos = eos_arr(0); disk_params.kbmu = constants.GetKBCode() / (mu * constants.GetAMUCode()); disk_params.pres_min = eos.PressureFromDensityInternalEnergy(disk_params.dens_min, disk_params.sie_min); @@ -410,7 +411,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - const auto &eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Disk parameters auto disk_params = artemis_pkg->Param("disk_params"); @@ -443,9 +444,8 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { pmb->par_for( "disk", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { - DiskICImpl(v, 0, k, j, i, pco, eos_d, dp, particles, npart); + DiskICImpl(v, 0, k, j, i, pco, eos_d(0), dp, particles, npart); }); - if (dp.do_imc) jaybenne::InitializeRadiation(md.get(), true); } //---------------------------------------------------------------------------------------- @@ -471,7 +471,7 @@ void DiskBoundaryVisc(std::shared_ptr> &mbd, bool coarse) { // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Packing static auto descriptors = ArtemisUtils::GetBoundaryPackDescriptorMap< @@ -566,8 +566,8 @@ void DiskBoundaryVisc(std::shared_ptr> &mbd, bool coarse) { const Real dx = std::log(xvp1[ix1] / xvm1[ix1]); const Real xmadx = xma / dx; - const Real nua = ViscosityProfile(dp, eos_d, xcyla[0], xcyla[2]); - const Real nug = ViscosityProfile(dp, eos_d, xcyl[0], xcyl[2]); + const Real nua = ViscosityProfile(dp, eos_d(0), xcyla[0], xcyla[2]); + const Real nug = ViscosityProfile(dp, eos_d(0), xcyl[0], xcyl[2]); // Viscous BC for gas if (do_gas) { @@ -690,7 +690,7 @@ void DiskBoundaryIC(std::shared_ptr> &mbd, bool coarse) { auto disk_params = artemis_pkg->Param("disk_params"); auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); static auto descriptors = ArtemisUtils::GetBoundaryPackDescriptorMap< gas::prim::density, gas::prim::velocity, gas::prim::sie, dust::prim::density, @@ -714,7 +714,7 @@ void DiskBoundaryIC(std::shared_ptr> &mbd, bool coarse) { pmb->par_for_bndry( "DiskInnerX1", nb, BDY, parthenon::TopologicalElement::CC, coarse, fine, KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { - DiskICImpl(v, 0, k, j, i, pco, eos_d, dp, particles, npart); + DiskICImpl(v, 0, k, j, i, pco, eos_d(0), dp, particles, npart); }); } @@ -736,7 +736,7 @@ void DiskBoundaryExtrap(std::shared_ptr> &mbd, bool coarse) // Extract gas parameters auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); auto disk_params = artemis_pkg->Param("disk_params"); auto &dp = disk_params; diff --git a/src/pgen/escape_1d.hpp b/src/pgen/escape_1d.hpp new file mode 100644 index 00000000..6fb57316 --- /dev/null +++ b/src/pgen/escape_1d.hpp @@ -0,0 +1,244 @@ +//======================================================================================== +// (C) (or copyright) 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +#ifndef PGEN_ESCAPE_1D_HPP_ +#define PGEN_ESCAPE_1D_HPP_ +//! \file escape_1d.hpp +//! \brief 1D spherical planetary atmospheric escape test. +//! +//! ## Problem Description +//! Two (or more) gas species occupy a 1D spherical atmosphere under the gravity of a +//! central point mass. Species 0 is the "light" component (e.g., hydrogen, mu=1) that +//! has enough thermal energy to escape via a Parker-type wind. Species 1 is the +//! "heavy" component (e.g., oxygen, mu=16) that is gravitationally bound. +//! +//! Chapman-Cowling drag couples the two fluids. The test demonstrates: +//! (a) Qualitative: heavy species is entrained and dragged outward by light species. +//! (b) Quantitative: mass flux of light species matches the isothermal Parker wind +//! solution at late times (in the absence of heavy species). +//! +//! ## Initial Conditions +//! Isothermal hydrostatic equilibrium initialised from the stellar surface at r0. +//! Each species satisfies: +//! rho_n(r) = rho0_n * exp(-mu_n * G*M / (R_gas_n * T0) * (1/r0 - 1/r)) +//! where R_gas_n = kB / (mu_n * m_amu). +//! +//! ## Boundary Conditions +//! - inner_x1: hydrostatic inflow (density fixed, velocity set by Parker wind) +//! - outer_x1: supersonic outflow (zero-gradient extrapolation) +//! +//! ## Required input parameters (in [problem] block) +//! - rho0_0 : base density species 0 [code units] +//! - rho0_1 : base density species 1 [code units] +//! - T0 : isothermal temperature [K] +//! - r0 : inner radius (base of atmosphere) [code units] +//! - gm : G*M [code units] (central body, can match gravity package) +//! - nspecies : number of species (2 for this test) + +// C/C++ headers +#include +#include +#include + +// Artemis headers +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/eos/eos.hpp" + +using ArtemisUtils::EOS; +using ArtemisUtils::VI; + +namespace escape_1d { + +struct Escape1DParams { + int nspecies; + Real rho0[8]; // base densities (up to 8 species) + Real T0; // common isothermal temperature [K] + Real r0; // inner boundary radius [code length] + Real gm; // G*M [code length^3 / code time^2] + Real kbmu[8]; // kB/(mu_n * m_amu) = R_gas per species [code velocity^2 / K] + Real mu[8]; // molecular masses [AMU] + Real dfloor; + Real siefloor; +}; + +//---------------------------------------------------------------------------------------- +//! \fn void InitEscape1DParams +//! \brief Reads problem parameters and stores in artemis package. +inline void InitEscape1DParams(MeshBlock *pmb, ParameterInput *pin) { + auto &artemis_pkg = pmb->packages.Get("artemis"); + Params ¶ms = artemis_pkg->AllParams(); + if (params.hasKey("escape1d_params")) return; + + auto &gas_pkg = pmb->packages.Get("gas"); + const Real kb = gas_pkg->Param("kb"); + const Real amu = gas_pkg->Param("amu"); + const int nspecies = gas_pkg->Param("nspecies"); + PARTHENON_REQUIRE(nspecies >= 2, "escape_1d requires at least 2 gas species"); + PARTHENON_REQUIRE(nspecies <= 8, "escape_1d supports at most 8 species"); + + // Per-species molecular masses from gas/species_mu + auto mu_arr = gas_pkg->Param>("mu"); + auto mu_h = mu_arr.GetHostMirrorAndCopy(); + + Escape1DParams ep{}; + ep.nspecies = nspecies; + ep.T0 = pin->GetOrAddReal("problem", "T0", 1.0e4); + ep.r0 = pin->GetOrAddReal("problem", "r0", 1.0); + ep.gm = pin->GetOrAddReal("problem", "gm", 1.0); + ep.dfloor = gas_pkg->Param("dfloor"); + ep.siefloor = gas_pkg->Param("siefloor"); + + const std::string rho0_key = "rho0_"; + for (int n = 0; n < nspecies; ++n) { + ep.mu[n] = mu_h(n); + ep.kbmu[n] = kb / (ep.mu[n] * amu); + ep.rho0[n] = pin->GetOrAddReal("problem", rho0_key + std::to_string(n), 1.0e-4); + } + params.Add("escape1d_params", ep); +} + +// --------------------------------------------------------------------------- +// Helpers: isothermal hydrostatic density at radius r +// --------------------------------------------------------------------------- +KOKKOS_FORCEINLINE_FUNCTION +Real IsothermalDens(const Escape1DParams &ep, const int n, const Real r) { + const Real H_inv = ep.gm / (ep.kbmu[n] * ep.T0); // inverse scale height at r0 + const Real exponent = -H_inv * (1.0 / ep.r0 - 1.0 / r); + return ep.rho0[n] * std::exp(exponent); +} + +//---------------------------------------------------------------------------------------- +//! \fn void ProblemGenerator +//! \brief Sets initial conditions for the atmospheric escape problem. +template +inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + + auto artemis_pkg = pmb->packages.Get("artemis"); + const bool do_gas = artemis_pkg->Param("do_gas"); + PARTHENON_REQUIRE(do_gas, "escape_1d requires gas hydrodynamics!"); + + auto gas_pkg = pmb->packages.Get("gas"); + const auto &eos_d = gas_pkg->Param>("eos_d"); + + // Ensure params are initialised (called here and also via InitMeshBlockData) + InitEscape1DParams(pmb, pin); + const auto ep = artemis_pkg->Param("escape1d_params"); + const int ns = ep.nspecies; + + auto &md = pmb->meshblock_data.Get(); + for (auto &var : md->GetVariableVector()) { + if (!var->IsAllocated()) pmb->AllocateSparse(var->label()); + } + + static auto desc = + MakePackDescriptor( + (pmb->resolved_packages).get()); + auto v = desc.GetPack(md.get()); + + static auto desc_g = MakePackDescriptor((pmb->resolved_packages).get()); + auto vg = desc_g.GetPack(md.get()); + + const auto &cpars = artemis_pkg->template Param("coord_params"); + auto &pco = pmb->coords; + + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + pmb->par_for( + "escape_1d_ic", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + geometry::Coords coords(cpars, pco, k, j, i); + const Real r = coords.x1v(); + for (int n = 0; n < ns; ++n) { + const Real rho = std::max(IsothermalDens(ep, n, r), ep.dfloor); + const Real sie = std::max( + eos_d(n).InternalEnergyFromDensityTemperature(rho, ep.T0), ep.siefloor); + v(0, gas::prim::density(n), k, j, i) = rho; + // Initialize at rest; the Parker wind will develop from the inner BC + for (int d = 0; d < 3; ++d) + v(0, gas::prim::velocity(VI(n, d)), k, j, i) = 0.0; + v(0, gas::prim::sie(n), k, j, i) = sie; + } + }); // par_for IC +} // ProblemGenerator + +//---------------------------------------------------------------------------------------- +//! \fn Escape1DInnerX1 +//! \brief Inner boundary: hydrostatic density fix, outward velocity extrapolation. +//! For the Parker wind test the inner boundary is below the sonic point so the +//! flow is subsonic and we fix density/temperature while allowing the velocity to +//! be set by the solution in the ghost cells (extrapolation inwards). +template +void Escape1DInnerX1(std::shared_ptr> &mbd, bool coarse) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + auto pmb = mbd->GetBlockPointer(); + auto artemis_pkg = pmb->packages.Get("artemis"); + const auto ep = artemis_pkg->Param("escape1d_params"); + const int ns = ep.nspecies; + + auto gas_pkg = pmb->packages.Get("gas"); + const auto eos_d = gas_pkg->Param>("eos_d"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + auto &pco = pmb->coords; + + static auto desc = + MakePackDescriptor((pmb->resolved_packages).get()); + auto v = desc.GetPack(mbd.get()); + + const IndexRange ib = pmb->cellbounds.GetBoundsI(BND, parthenon::TE::CC); + const IndexRange jb = pmb->cellbounds.GetBoundsJ(BND, parthenon::TE::CC); + const IndexRange kb = pmb->cellbounds.GetBoundsK(BND, parthenon::TE::CC); + const IndexRange ibi = + pmb->cellbounds.GetBoundsI(parthenon::IndexDomain::interior, parthenon::TE::CC); + + pmb->par_for( + "escape_1d_innerX1", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // Mirror index into interior + const int im = ibi.s + (ibi.s - i - 1); + const int isafe = std::max(im, ibi.s); + + geometry::Coords coords(cpars, pco, k, j, i); + const Real r = coords.x1v(); + + for (int n = 0; n < ns; ++n) { + const Real rho = std::max(IsothermalDens(ep, n, r), ep.dfloor); + const Real sie = std::max( + eos_d(n).InternalEnergyFromDensityTemperature(rho, ep.T0), ep.siefloor); + + // Velocity: extrapolate from interior (allow Parker wind to develop) + const Real rho_int = v(0, gas::cons::density(n), k, j, isafe); + const Real vr_int = + (rho_int > 0.0) ? v(0, gas::cons::momentum(VI(n, 0)), k, j, isafe) / rho_int + : 0.0; + const Real vr = std::max(vr_int, 0.0); // only outward flow at inner BC + + v(0, gas::cons::density(n), k, j, i) = rho; + v(0, gas::cons::momentum(VI(n, 0)), k, j, i) = rho * vr; + v(0, gas::cons::momentum(VI(n, 1)), k, j, i) = 0.0; + v(0, gas::cons::momentum(VI(n, 2)), k, j, i) = 0.0; + const Real ke = 0.5 * rho * vr * vr; + v(0, gas::cons::internal_energy(n), k, j, i) = rho * sie; + v(0, gas::cons::total_energy(n), k, j, i) = rho * sie + ke; + } + }); +} +} // namespace escape_1d + +#endif // PGEN_ESCAPE_1D_HPP_ diff --git a/src/pgen/gaussian_bump.hpp b/src/pgen/gaussian_bump.hpp index 8ff99c27..9744f739 100644 --- a/src/pgen/gaussian_bump.hpp +++ b/src/pgen/gaussian_bump.hpp @@ -65,7 +65,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { bump_params.vfac = pin->GetOrAddReal("problem", "vx2_bump", 0.0); bump_params.wfac = pin->GetOrAddReal("problem", "vx3_bump", 0.0); - ArtemisUtils::EOS eos; + ParArray1D eos; if (do_gas) { auto gas_pkg = pmb->packages.Get("gas"); PARTHENON_REQUIRE((gas_pkg->Param("nspecies") == 1), @@ -75,7 +75,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { bump_params.g_vx2 = pin->GetOrAddReal("problem", "gas_vx2", 0.0); bump_params.g_vx3 = pin->GetOrAddReal("problem", "gas_vx3", 0.0); bump_params.g_pres = pin->GetOrAddReal("problem", "gas_pres", 1.0); - eos = gas_pkg->Param("eos_d"); + eos = gas_pkg->template Param>("eos_d"); } if (do_dust) { auto dust_pkg = pmb->packages.Get("dust"); @@ -179,16 +179,17 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (pars.tfac > 0.0) { // P = const, T = T0*(1 + f) // Input pressure and density to get background T - temp = ArtemisUtils::TofPR(eos, pres, dens); + temp = ArtemisUtils::TofPR(eos(0), pres, dens); temp *= (1. + pars.tfac * bump); // Input pressure and temperature - v(0, gas::prim::density(0), k, j, i) = ArtemisUtils::RofPT(eos, pres, temp); - v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPT(eos, pres, temp); + v(0, gas::prim::density(0), k, j, i) = + ArtemisUtils::RofPT(eos(0), pres, temp); + v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPT(eos(0), pres, temp); } else { dens = pars.g_rho * (1. + pars.dfac * bump); v(0, gas::prim::density(0), k, j, i) = dens; // input density and pressure - v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos, pres, dens); + v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos(0), pres, dens); } } if (do_dust) { diff --git a/src/pgen/kh.hpp b/src/pgen/kh.hpp index 28a6e4c5..2e805fa2 100644 --- a/src/pgen/kh.hpp +++ b/src/pgen/kh.hpp @@ -47,7 +47,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // Extract parameters from packages auto artemis_pkg = pmb->packages.Get("artemis"); auto gas_pkg = pmb->packages.Get("gas"); - const auto &eos = gas_pkg->Param("eos_d"); + const auto &eos = gas_pkg->template Param>("eos_d"); const bool do_gas = artemis_pkg->Param("do_gas"); KH_params.y1 = pin->GetOrAddReal("problem", "y1", 0.5); @@ -101,7 +101,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { std::exp(-SQR((zc - pars.y2) / pars.sigma))); v(0, gas::prim::density(0), k, j, i) = dens; - v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos, pars.pres0, dens); + v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos(0), pars.pres0, dens); v(0, gas::prim::velocity(0), k, j, i) = vx; v(0, gas::prim::velocity(1), k, j, i) = (!three_d) * vz; v(0, gas::prim::velocity(2), k, j, i) = three_d * vz; diff --git a/src/pgen/lw.hpp b/src/pgen/lw.hpp index 508fbaa3..d262af87 100644 --- a/src/pgen/lw.hpp +++ b/src/pgen/lw.hpp @@ -52,7 +52,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { lw_params.rho1 = pin->GetOrAddReal("problem", "rho1", 1.0); lw_params.pres1 = pin->GetOrAddReal("problem", "pres1", 1.0); - const auto &eos = gas_pkg->Param("eos_d"); + const auto &eos = gas_pkg->Param>("eos_d"); // packing and capture variables for kernel auto &md = pmb->meshblock_data.Get(); @@ -89,7 +89,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real pres = vf1 * pars.pres1 + vf2 * pars.pres0; v(0, gas::prim::density(0), k, j, i) = dens; - v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos, pres, dens); + v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos(0), pres, dens); v(0, gas::prim::velocity(0), k, j, i) = 0.0; v(0, gas::prim::velocity(1), k, j, i) = 0.0; v(0, gas::prim::velocity(2), k, j, i) = 0.0; diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 9ae84c8f..6cfe5cf8 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -26,6 +26,7 @@ #include "constant.hpp" #include "crooked_pipe.hpp" #include "disk.hpp" +#include "escape_1d.hpp" #include "gaussian_bump.hpp" #include "geometry/geometry.hpp" #include "grav_slab.hpp" @@ -78,6 +79,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { rt::ProblemGenerator(pmb, pin); } else if (name == "shock") { shock::ProblemGenerator(pmb, pin); + } else if (name == "escape_1d") { + escape_1d::ProblemGenerator(pmb, pin); } else if (name == "strat") { strat::ProblemGenerator(pmb, pin); } else if (name == "thermalization") { @@ -107,6 +110,8 @@ void InitMeshBlockData(MeshBlock *pmb, ParameterInput *pin) { disk::InitDiskParams(pmb, pin); } else if (name == "shock") { shock::InitShockParams(pmb, pin); + } else if (name == "escape_1d") { + escape_1d::InitEscape1DParams(pmb, pin); } else if (name == "strat") { strat::InitStratParams(pmb, pin); } diff --git a/src/pgen/problem_modifier.hpp b/src/pgen/problem_modifier.hpp index 162e3c02..2bd23287 100644 --- a/src/pgen/problem_modifier.hpp +++ b/src/pgen/problem_modifier.hpp @@ -122,6 +122,9 @@ void ProblemModifier(parthenon::ParthenonManager *pman) { strat::ExtrapInnerX3); pman->app_input->RegisterBoundaryCondition(BF::outer_x3, "extrap", strat::ExtrapOuterX3); + } else if (artemis_problem == "escape_1d") { + pman->app_input->RegisterBoundaryCondition( + BF::inner_x1, "hydrostatic", escape_1d::Escape1DInnerX1); } else if (artemis_problem == "beam") { pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "beam", beam::BeamInnerX1); diff --git a/src/pgen/rt.hpp b/src/pgen/rt.hpp index f2d6a8db..92b91d28 100644 --- a/src/pgen/rt.hpp +++ b/src/pgen/rt.hpp @@ -20,6 +20,7 @@ #include "utils/eos/eos.hpp" using ArtemisUtils::EOS; +using ArtemisUtils::VI; namespace { struct RTParams { @@ -57,7 +58,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { RT_params.freq = pin->GetOrAddReal("problem", "frequency", 6 * M_PI); RT_params.amp = pin->GetOrAddReal("problem", "amplitude", 0.01); - const auto &eos = gas_pkg->Param("eos_d"); + const auto &eos = gas_pkg->template Param>("eos_d"); // packing and capture variables for kernel auto &md = pmb->meshblock_data.Get(); @@ -111,11 +112,14 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real p0 = pars.pres0 + (upper)*gx * (pars.y0 - zmin) * pars.rho0; const Real pres = p0 + gx * (zc - z0) * dens; - v(0, gas::prim::density(0), k, j, i) = dens; - v(0, gas::prim::sie(0), k, j, i) = ArtemisUtils::EofPR(eos, pres, dens); - v(0, gas::prim::velocity(0), k, j, i) = 0.0; - v(0, gas::prim::velocity(1), k, j, i) = pars.amp * std::cos(pars.freq * xc); - v(0, gas::prim::velocity(2), k, j, i) = 0.0; + const int fid = (upper) ? 1 : 0; + + v(0, gas::prim::density(fid), k, j, i) = dens; + v(0, gas::prim::sie(fid), k, j, i) = ArtemisUtils::EofPR(eos(0), pres, dens); + v(0, gas::prim::velocity(VI(fid, 0)), k, j, i) = 0.0; + v(0, gas::prim::velocity(VI(fid, 1)), k, j, i) = + pars.amp * std::cos(pars.freq * xc); + v(0, gas::prim::velocity(VI(fid, 2)), k, j, i) = 0.0; }); } diff --git a/src/pgen/shock.hpp b/src/pgen/shock.hpp index 9be64caf..5354a07a 100644 --- a/src/pgen/shock.hpp +++ b/src/pgen/shock.hpp @@ -80,14 +80,11 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const bool do_moment = artemis_pkg->Param("do_moment"); PARTHENON_REQUIRE(do_gas, "The shock problem requires gas hydrodynamics!"); PARTHENON_REQUIRE(!(do_dust), "The shock problem does not permit dust hydrodynamics!"); - auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); - + const auto &eos_d = pmb->packages.Get("gas")->Param>("eos_d"); Real ar = Null(); if (do_moment) { ar = pmb->packages.Get("moments")->Param("arad"); } - - // packing and capture variables for kernel auto &md = pmb->meshblock_data.Get(); for (auto &var : md->GetVariableVector()) { if (!var->IsAllocated()) pmb->AllocateSparse(var->label()); @@ -123,7 +120,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(0, gas::prim::velocity(1), k, j, i) = 0.0; v(0, gas::prim::velocity(2), k, j, i) = 0.0; v(0, gas::prim::sie(0), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(rho, T); + eos_d(0).InternalEnergyFromDensityTemperature(rho, T); if (do_moment) { v(0, rad::prim::energy(0), k, j, i) = ar * SQR(SQR(T)); v(0, rad::prim::flux(0), k, j, i) = 0.0; @@ -147,7 +144,8 @@ inline void ShockInnerX1(std::shared_ptr> &mbd, bool coarse) auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_moment = artemis_pkg->Param("do_moment"); auto shkp = artemis_pkg->Param("shock_params"); - auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); + const auto &eos_d = pmb->packages.Get("gas")->Param>("eos_d"); + Real ar = Null(); if (do_moment) { ar = pmb->packages.Get("moments")->Param("arad"); @@ -171,7 +169,7 @@ inline void ShockInnerX1(std::shared_ptr> &mbd, bool coarse) v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = 0.0; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = 0.0; v(0, gas::prim::sie(n), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(shkp.rhol, shkp.tl); + eos_d(0).InternalEnergyFromDensityTemperature(shkp.rhol, shkp.tl); } if (do_moment) { for (int n = 0; n < v.GetSize(0, rad::prim::energy()); ++n) { @@ -198,7 +196,7 @@ inline void ShockOuterX1(std::shared_ptr> &mbd, bool coarse) auto artemis_pkg = pmb->packages.Get("artemis"); const bool do_moment = artemis_pkg->Param("do_moment"); auto shkp = artemis_pkg->Param("shock_params"); - auto eos_d = pmb->packages.Get("gas")->Param("eos_d"); + const auto &eos_d = pmb->packages.Get("gas")->Param>("eos_d"); Real ar = Null(); if (do_moment) { ar = pmb->packages.Get("moments")->Param("arad"); @@ -222,7 +220,7 @@ inline void ShockOuterX1(std::shared_ptr> &mbd, bool coarse) v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = 0.0; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = 0.0; v(0, gas::prim::sie(n), k, j, i) = - eos_d.InternalEnergyFromDensityTemperature(shkp.rhor, shkp.tr); + eos_d(0).InternalEnergyFromDensityTemperature(shkp.rhor, shkp.tr); } if (do_moment) { for (int n = 0; n < v.GetSize(0, rad::prim::energy()); ++n) { diff --git a/src/pgen/strat.hpp b/src/pgen/strat.hpp index 6e5f4a3a..c7631206 100644 --- a/src/pgen/strat.hpp +++ b/src/pgen/strat.hpp @@ -86,7 +86,8 @@ inline void InitStratParams(MeshBlock *pmb, ParameterInput *pin) { auto &gas_pkg = pmb->packages.Get("gas"); const auto mu = gas_pkg->Param("mu"); - const auto eos = gas_pkg->Param("eos_h"); + const auto eos_arr = gas_pkg->Param>("eos_h"); + const EOS &eos = eos_arr(0); auto &constants = artemis_pkg->Param("constants"); strat_params.kbmu = constants.GetKBCode() / (mu * constants.GetAMUCode()); strat_params.dfloor = gas_pkg->Param("dfloor"); @@ -152,7 +153,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto strat_params = artemis_pkg->Param("strat_params"); auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); const bool is_ideal = (gas_pkg->template Param("eos_type") == "ideal"); // Dimensionality @@ -189,9 +190,9 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real vx3 = 0.0; const Real temp = pars.temp0; const Real dens = - is_ideal ? InitialDensity(pars, z) : InitialDensity(eos_d, pars, z); - const Real sie = std::max(pars.siefloor, - eos_d.InternalEnergyFromDensityTemperature(dens, temp)); + is_ideal ? InitialDensity(pars, z) : InitialDensity(eos_d(0), pars, z); + const Real sie = std::max( + pars.siefloor, eos_d(0).InternalEnergyFromDensityTemperature(dens, temp)); v(0, gas::prim::density(0), k, j, i) = dens; v(0, gas::prim::velocity(0), k, j, i) = vx1; @@ -445,7 +446,7 @@ inline void ShearInnerX2(std::shared_ptr> &mbd, bool coarse) // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Coordinates and indexing const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; @@ -479,10 +480,11 @@ inline void ShearInnerX2(std::shared_ptr> &mbd, bool coarse) const Real vx3g = outflow ? gv3 : 0.0; const Real densg = outflow ? v(0, gas::prim::density(n), k, js, i) : InitialDensity(pars, z); - const Real sieg = outflow ? v(0, gas::prim::sie(n), k, js, i) - : std::max(pars.siefloor, - eos_d.InternalEnergyFromDensityTemperature( - densg, pars.temp0)); + const Real sieg = + outflow + ? v(0, gas::prim::sie(n), k, js, i) + : std::max(pars.siefloor, eos_d(n).InternalEnergyFromDensityTemperature( + densg, pars.temp0)); v(0, gas::prim::velocity(VI(n, 0)), k, j, i) = vx1g; v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = vx2g; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = vx3g; @@ -566,7 +568,7 @@ inline void ShearOuterX2(std::shared_ptr> &mbd, bool coarse) // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Coordinates and indexing const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; @@ -600,10 +602,11 @@ inline void ShearOuterX2(std::shared_ptr> &mbd, bool coarse) const Real vx3g = outflow ? gv3 : 0.0; const Real densg = outflow ? v(0, gas::prim::density(n), k, je, i) : InitialDensity(pars, z); - const Real sieg = outflow ? v(0, gas::prim::sie(n), k, je, i) - : std::max(pars.siefloor, - eos_d.InternalEnergyFromDensityTemperature( - densg, pars.temp0)); + const Real sieg = + outflow + ? v(0, gas::prim::sie(n), k, je, i) + : std::max(pars.siefloor, eos_d(n).InternalEnergyFromDensityTemperature( + densg, pars.temp0)); v(0, gas::prim::velocity(VI(n, 0)), k, j, i) = vx1g; v(0, gas::prim::velocity(VI(n, 1)), k, j, i) = vx2g; v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = vx3g; @@ -674,7 +677,7 @@ inline void ExtrapInnerX3(std::shared_ptr> &mbd, bool coarse // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Coordinates and indexing const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; @@ -705,9 +708,9 @@ inline void ExtrapInnerX3(std::shared_ptr> &mbd, bool coarse const Real vx3g = (gv3 > 0.0) ? 0.0 : gv3; const Real &gd = v(0, gas::prim::density(n), ks, j, i); const Real &gsie = v(0, gas::prim::sie(n), ks, j, i); - const Real Tg = eos_d.TemperatureFromDensityInternalEnergy(gd, gsie); - const Real pm = eos_d.PressureFromDensityTemperature(gd * (1. - 1e-6), Tg); - const Real pp = eos_d.PressureFromDensityTemperature(gd * (1. + 1e-6), Tg); + const Real Tg = eos_d(n).TemperatureFromDensityInternalEnergy(gd, gsie); + const Real pm = eos_d(n).PressureFromDensityTemperature(gd * (1. - 1e-6), Tg); + const Real pp = eos_d(n).PressureFromDensityTemperature(gd * (1. + 1e-6), Tg); const Real dPdrho = (pp - pm) / (gd * 1e-6); const Real efac = std::exp(-(SQR(z) - SQR(z0)) * SQR(pars.Om0) / (2.0 * dPdrho + Fuzz())); @@ -718,7 +721,7 @@ inline void ExtrapInnerX3(std::shared_ptr> &mbd, bool coarse v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = vx3g; v(0, gas::prim::density(n), k, j, i) = rhog; v(0, gas::prim::sie(n), k, j, i) = std::max( - pars.siefloor, eos_d.InternalEnergyFromDensityTemperature(rhog, Tg)); + pars.siefloor, eos_d(n).InternalEnergyFromDensityTemperature(rhog, Tg)); } if (pars.do_dust) { @@ -776,7 +779,7 @@ inline void ExtrapOuterX3(std::shared_ptr> &mbd, bool coarse // Extract gas package and params auto &gas_pkg = pmb->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); // Coordinates and indexing const auto &pco = (coarse) ? pmb->pmr->GetCoarseCoords() : pmb->coords; @@ -807,9 +810,9 @@ inline void ExtrapOuterX3(std::shared_ptr> &mbd, bool coarse const Real vx3g = (gv3 < 0.0) ? 0.0 : gv3; const Real &gd = v(0, gas::prim::density(n), ke, j, i); const Real &gsie = v(0, gas::prim::sie(n), ke, j, i); - const Real Tg = eos_d.TemperatureFromDensityInternalEnergy(gd, gsie); - const Real pm = eos_d.PressureFromDensityTemperature(gd * (1. - 1e-6), Tg); - const Real pp = eos_d.PressureFromDensityTemperature(gd * (1. + 1e-6), Tg); + const Real Tg = eos_d(n).TemperatureFromDensityInternalEnergy(gd, gsie); + const Real pm = eos_d(n).PressureFromDensityTemperature(gd * (1. - 1e-6), Tg); + const Real pp = eos_d(n).PressureFromDensityTemperature(gd * (1. + 1e-6), Tg); const Real dPdrho = (pp - pm) / (gd * 1e-6); const Real efac = std::exp(-(SQR(z) - SQR(z0)) * SQR(pars.Om0) / (2.0 * dPdrho + Fuzz())); @@ -820,7 +823,7 @@ inline void ExtrapOuterX3(std::shared_ptr> &mbd, bool coarse v(0, gas::prim::velocity(VI(n, 2)), k, j, i) = vx3g; v(0, gas::prim::density(n), k, j, i) = rhog; v(0, gas::prim::sie(n), k, j, i) = std::max( - pars.siefloor, eos_d.InternalEnergyFromDensityTemperature(rhog, Tg)); + pars.siefloor, eos_d(n).InternalEnergyFromDensityTemperature(rhog, Tg)); } if (pars.do_dust) { diff --git a/src/pgen/thermalization.hpp b/src/pgen/thermalization.hpp index 2079d377..dad2037a 100644 --- a/src/pgen/thermalization.hpp +++ b/src/pgen/thermalization.hpp @@ -47,7 +47,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(do_gas, "Thermalization problem requires gas!"); PARTHENON_REQUIRE(!(do_dust), "Thermalization problem does not permit dust!"); auto gas_pkg = pmb->packages.Get("gas"); - const auto eos = gas_pkg->Param("eos_d"); + const auto eos = gas_pkg->Param>("eos_d"); Real ar = Null(); if (do_moment) { auto rad_pkg = pmb->packages.Get("moments"); @@ -83,7 +83,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { KOKKOS_LAMBDA(const int k, const int j, const int i) { v(0, gas::prim::density(), k, j, i) = rho; v(0, gas::prim::sie(), k, j, i) = - eos.InternalEnergyFromDensityTemperature(rho, trad); + eos(0).InternalEnergyFromDensityTemperature(rho, trad); }); jaybenne::InitializeRadiation(md.get(), true); } @@ -97,7 +97,7 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(0, gas::prim::velocity(1), k, j, i) = 0.0; v(0, gas::prim::velocity(2), k, j, i) = 0.0; v(0, gas::prim::sie(), k, j, i) = - eos.InternalEnergyFromDensityTemperature(rho, tgas); + eos(0).InternalEnergyFromDensityTemperature(rho, tgas); if (do_moment) { v(0, rad::prim::energy(), k, j, i) = ar * SQR(SQR(trad)); diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index c471f3a4..adbf7cb1 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -41,7 +41,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); auto opac_d = gas_pkg->template Param("opacity_d"); auto scat_d = gas_pkg->template Param("scattering_d"); auto dflr = gas_pkg->template Param("dfloor"); @@ -118,16 +118,17 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { // Note(AMD): There is some floating point difference between the internal energy // used to compute the temperature and the internal energy obtained from that - // temperature: T = eos_d.TemperatureFromDensityInternalEnergy(dens, eg/dens); eg - // /= dens * eos_d.InternalEnergyFromDensityTemperature(dens,T) + // temperature: T = eos_d(0).TemperatureFromDensityInternalEnergy(dens, eg/dens); + // eg + // /= dens * eos_d(0).InternalEnergyFromDensityTemperature(dens,T) // // Because of this, zero opacity problems will not result in zero change as // expected. Thus, we recalculate the internal and total energies from the // temperature. This does not affect energy conservation because at the end of the // step we update the energy with an increment. - Real T = eos_d.TemperatureFromDensityInternalEnergy(dens, e0 / dens); - e0 = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T); + Real T = eos_d(0).TemperatureFromDensityInternalEnergy(dens, e0 / dens); + e0 = dens * eos_d(0).InternalEnergyFromDensityTemperature(dens, T); Real e = e0; Real B = arad * SQR(SQR(T)); @@ -140,8 +141,8 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { Real inner_err = 0.; for (inner_iter = 0; inner_iter < inner_max; inner_iter++) { T = std::pow(B / arad, 0.25); - e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; - const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); + e = eos_d(0).InternalEnergyFromDensityTemperature(dens, T) * dens; + const Real Cv = dens * eos_d(0).SpecificHeatFromDensityTemperature(dens, T); const Real a = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); const Real fleck = FleckFactor(arad, T, Cv); @@ -168,7 +169,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { PARTHENON_FAIL("Radiation matter coupling did not converge!"); } T = std::pow(B / arad, 0.25); - e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; + e = eos_d(0).InternalEnergyFromDensityTemperature(dens, T) * dens; const Real dEg = e - e0; Real a = chat * dt * (opac_d.RosselandMeanAbsorptionCoefficient(dens, T) + @@ -210,7 +211,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); auto opac_d = gas_pkg->template Param("opacity_d"); auto scat_d = gas_pkg->template Param("scattering_d"); auto dflr = gas_pkg->template Param("dfloor"); @@ -275,16 +276,17 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Note(AMD): There is some floating point difference between the internal energy // used to compute the temperature and the internal energy obtained from that - // temperature: T = eos_d.TemperatureFromDensityInternalEnergy(dens, eg/dens); eg - // /= dens * eos_d.InternalEnergyFromDensityTemperature(dens,T) + // temperature: T = eos_d(0).TemperatureFromDensityInternalEnergy(dens, eg/dens); + // eg + // /= dens * eos_d(0).InternalEnergyFromDensityTemperature(dens,T) // // Because of this, zero opacity problems will not result in zero change as // expected. Thus, we recalculate the internal and total energies from the // temperature. This does not affect energy conservation because at the end of the // step we update the energy with an increment. - Real T = eos_d.TemperatureFromDensityInternalEnergy( + Real T = eos_d(0).TemperatureFromDensityInternalEnergy( dens, v0(b, gas::cons::internal_energy(), k, j, i) / dens); - Real eg0 = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T); + Real eg0 = dens * eos_d(0).InternalEnergyFromDensityTemperature(dens, T); Real B = arad * SQR(SQR(T)); const auto vb = RotatingFrame::BackgroundVelocity( @@ -359,9 +361,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // start inner iteration for (B,E) for (inner_iter = 1; inner_iter <= inner_max; inner_iter++) { T = std::pow(eref * B / arad, 0.25); - Real eint = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; + Real eint = + dens * eos_d(0).InternalEnergyFromDensityTemperature(dens, T) / eref; Real et = ke + eint; - const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); + const Real Cv = dens * eos_d(0).SpecificHeatFromDensityTemperature(dens, T); const Real fleck = FleckFactor(arad, T, Cv); const Real sigp = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); @@ -404,7 +407,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Have new E and T T = std::pow(eref * B / arad, 0.25); - Real eg = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; + Real eg = dens * eos_d(0).InternalEnergyFromDensityTemperature(dens, T) / eref; dEg = eg - eg0; const Real sigp = diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 9a10d34d..61730f9a 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -90,9 +90,9 @@ TaskStatus SetOpacities(MeshData *md) { auto &resolved_pkgs = pm->resolved_packages; auto &gas_pkg = pm->packages.Get("gas"); - EOS eos_d = gas_pkg->template Param("eos_d"); - MeanOpacity opacity_d = gas_pkg->template Param("opacity_d"); - MeanScattering scattering_d = gas_pkg->template Param("scattering_d"); + const auto &eos_d = gas_pkg->template Param>("eos_d"); + auto &opacity_d = gas_pkg->template Param("opacity_d"); + auto &scattering_d = gas_pkg->template Param("scattering_d"); // Packing and indexing // TODO(): Will eventually incorporate other fluids @@ -111,7 +111,7 @@ TaskStatus SetOpacities(MeshData *md) { KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { const Real &rho = vmesh(b, gas::prim::density(), k, j, i); const Real &sie = vmesh(b, gas::prim::sie(), k, j, i); - const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); + const Real temp = eos_d(0).TemperatureFromDensityInternalEnergy(rho, sie); Real &aa = vmesh(b, rad::opac::absorption(), k, j, i); Real &ss = vmesh(b, rad::opac::scattering(), k, j, i); diff --git a/src/radiation/raytrace/raytrace.cpp b/src/radiation/raytrace/raytrace.cpp index 3f888bba..96412821 100644 --- a/src/radiation/raytrace/raytrace.cpp +++ b/src/radiation/raytrace/raytrace.cpp @@ -275,9 +275,9 @@ TaskStatus EvalOpac(MeshData *md) { auto &rt_pkg = pm->packages.Get("raytrace"); auto &gas_pkg = pm->packages.Get("gas"); - ArtemisUtils::EOS eos_d = gas_pkg->template Param("eos_d"); - ArtemisUtils::Opacity opacity_d = - rt_pkg->template Param("opacity_d"); + ParArray1D eos_d = + gas_pkg->Param>("eos_d"); + ArtemisUtils::Opacity opacity_d = rt_pkg->Param("opacity_d"); static auto desc = MakePackDescriptor( @@ -297,7 +297,7 @@ TaskStatus EvalOpac(MeshData *md) { //%%%%%%%%%%%%%%%% // Evaluated at T* //%%%%%%%%%%%%%%%% - const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); + const Real temp = eos_d(0).TemperatureFromDensityInternalEnergy(rho, sie); vmesh(b, rad::star::absorption(), k, j, i) = opacity_d.AbsorptionCoefficient(rho, temp, 1.0); }); diff --git a/src/utils/artemis_utils.cpp b/src/utils/artemis_utils.cpp index 5b014c38..e6c4b290 100644 --- a/src/utils/artemis_utils.cpp +++ b/src/utils/artemis_utils.cpp @@ -51,18 +51,19 @@ void PrintArtemisConfiguration(Packages_t &packages) { if (params.Get("do_moment")) msg += hfill + "Moment radiation\n"; printf("\n=====================================================\n"); printf(" ARTEMIS\n"); - printf(" name: %s\n", params.Get("job_name").c_str()); - printf(" problem: %s\n", params.Get("pgen_name").c_str()); - printf(" coordinates: %dD %s\n", nd, params.Get("coord_sys").c_str()); - printf(" integrator: %s\n", params.Get("integrator").c_str()); - printf(" MPI ranks: %d\n", parthenon::Globals::nranks); - printf(" dimensions: %dx%dx%d\n", nx[0], nx[1], nx[2]); - printf(" meshblock: %dx%dx%d\n", nb[0], nb[1], nb[2]); - printf(" Unit System: %s\n", units.GetSystemName().c_str()); - printf(" [L] = %.2e\n", units.GetLengthCodeToPhysical()); - printf(" [M] = %.2e\n", units.GetMassCodeToPhysical()); - printf(" [T] = %.2e\n", units.GetTimeCodeToPhysical()); - printf(" [K] = %.2e\n", units.GetTemperatureCodeToPhysical()); + printf(" name: %s\n", params.Get("job_name").c_str()); + printf(" problem: %s\n", params.Get("pgen_name").c_str()); + printf(" sparse: %s\n", params.Get("do_sparse") ? "yes" : "no"); + printf(" coordinates: %dD %s\n", nd, params.Get("coord_sys").c_str()); + printf(" integrator: %s\n", params.Get("integrator").c_str()); + printf(" MPI ranks: %d\n", parthenon::Globals::nranks); + printf(" dimensions: %dx%dx%d\n", nx[0], nx[1], nx[2]); + printf(" meshblock: %dx%dx%d\n", nb[0], nb[1], nb[2]); + printf(" Unit System: %s\n", units.GetSystemName().c_str()); + printf(" [L] = %.2e\n", units.GetLengthCodeToPhysical()); + printf(" [M] = %.2e\n", units.GetMassCodeToPhysical()); + printf(" [T] = %.2e\n", units.GetTimeCodeToPhysical()); + printf(" [K] = %.2e\n", units.GetTemperatureCodeToPhysical()); printf(" Active physics: %s", msg.c_str()); if (params.Get("do_nbody")) { auto nbody_pkg = packages.Get("nbody"); diff --git a/src/utils/diffusion/diffusion.hpp b/src/utils/diffusion/diffusion.hpp index f5273ec7..e5414f11 100644 --- a/src/utils/diffusion/diffusion.hpp +++ b/src/utils/diffusion/diffusion.hpp @@ -63,8 +63,8 @@ TaskStatus ZeroDiffusionImpl(MeshData *md, SparsePackFlux vf) { //! \brief Computes diffusion limited timestep template -Real EstimateTimestep(MeshData *md, DiffCoeffParams &dp, PKG &pkg, const EOS &eos, - SparsePackPrim vprim) { +Real EstimateTimestep(MeshData *md, DiffCoeffParams &dp, PKG &pkg, + const ParArray1D &eos, SparsePackPrim vprim) { PARTHENON_INSTRUMENT using parthenon::MakePackDescriptor; auto pm = md->GetParentPointer(); @@ -99,9 +99,9 @@ Real EstimateTimestep(MeshData *md, DiffCoeffParams &dp, PKG &pkg, const E // Get the maximum diffusion coefficient (if there's more than one) DiffusionCoeff diffcoeff; - Real mu = diffcoeff.Get(dp, coords, xv, dens, sie, eos); + Real mu = diffcoeff.Get(dp, coords, xv, dens, sie, eos(n)); if constexpr (DIFF == DiffType::conductivity_plaw) { - mu /= (dens * eos.SpecificHeatFromDensityInternalEnergy(dens, sie)); + mu /= (dens * eos(n).SpecificHeatFromDensityInternalEnergy(dens, sie)); } else if constexpr ((DIFF == DiffType::viscosity_plaw) || (DIFF == DiffType::viscosity_alpha)) { mu *= (1.0 + (dp.eta > 1.0) * (dp.eta - 1.0)) / dens; diff --git a/src/utils/diffusion/momentum_diffusion.hpp b/src/utils/diffusion/momentum_diffusion.hpp index 3fe60b7e..a6e76bce 100644 --- a/src/utils/diffusion/momentum_diffusion.hpp +++ b/src/utils/diffusion/momentum_diffusion.hpp @@ -600,7 +600,7 @@ TaskStatus MomentumFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, "Momentum diffusion only works with a gas fluid"); auto pm = md->GetParentPointer(); - auto eos_d = pkg->template Param("eos_d"); + const auto &eos_d = pkg->template Param>("eos_d"); Real qshear = 0.0, om0 = 0.0, gm_bg = 0.0; const bool do_oa = @@ -660,7 +660,7 @@ TaskStatus MomentumFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, multi_d, three_d, vprim, vg, divu); // 3. Viscosity values. No barrier DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il - 1, iu, vprim, eos_d, mu); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il - 1, iu, vprim, eos_d(n), mu); mbr.team_barrier(); @@ -713,7 +713,7 @@ TaskStatus MomentumFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, // 3. Viscosity values. No barrier DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d, mu_jm1); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d(n), mu_jm1); mbr.team_barrier(); if (j > jl) { @@ -767,7 +767,7 @@ TaskStatus MomentumFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, // 2. Viscosity values. No barrier DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d, mu_km1); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d(n), mu_km1); mbr.team_barrier(); if (k > kl) { diff --git a/src/utils/diffusion/thermal_diffusion.hpp b/src/utils/diffusion/thermal_diffusion.hpp index 12d75014..a07bcb67 100644 --- a/src/utils/diffusion/thermal_diffusion.hpp +++ b/src/utils/diffusion/thermal_diffusion.hpp @@ -44,7 +44,7 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, "thermal diffusion only works with a gas fluid"); auto pm = md->GetParentPointer(); - auto eos_d = pkg->template Param("eos_d"); + const auto &eos_d = pkg->template Param>("eos_d"); static auto desc_g = MakePackDescriptor((pm->resolved_packages).get()); @@ -81,7 +81,7 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, // this returns conductivity not thermal diffusivity DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il - 1, iu, vprim, eos_d, kappa); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il - 1, iu, vprim, eos_d(n), kappa); mbr.team_barrier(); @@ -93,11 +93,11 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, const auto &xv_m = coords.GetCellCenter(vg, b, k, j, i - 1); const Real dx1 = coords.Distance(xv, xv_m); - const Real T = eos_d.TemperatureFromDensityInternalEnergy( + const Real T = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k, j, i), vprim(b, gas::prim::sie(n), k, j, i)); - const Real Tm = eos_d.TemperatureFromDensityInternalEnergy( + const Real Tm = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k, j, i - 1), vprim(b, gas::prim::sie(n), k, j, i - 1)); @@ -135,7 +135,7 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, } DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d, kappa); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d(n), kappa); mbr.team_barrier(); if (j > jl) { @@ -147,11 +147,11 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, const auto &xv_m = coords.GetCellCenter(vg, b, k, j - 1, i); const Real dx2 = coords.Distance(xv, xv_m); - const Real T = eos_d.TemperatureFromDensityInternalEnergy( + const Real T = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k, j, i), vprim(b, gas::prim::sie(n), k, j, i)); - const Real Tm = eos_d.TemperatureFromDensityInternalEnergy( + const Real Tm = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k, j - 1, i), vprim(b, gas::prim::sie(n), k, j - 1, i)); @@ -194,7 +194,7 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, // 2. Viscosity values. No barrier DiffusionCoeff diffcoeff; - diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d, kappa_km1); + diffcoeff.evaluate(dp, mbr, b, n, k, j, il, iu, vprim, eos_d(n), kappa_km1); mbr.team_barrier(); if (k > kl) { @@ -206,11 +206,11 @@ TaskStatus ThermalFluxImpl(MeshData *md, DiffCoeffParams dp, PKG &pkg, const auto &xv_m = coords.GetCellCenter(vg, b, k - 1, j, i); const Real dx3 = coords.Distance(xv, xv_m); - const Real T = eos_d.TemperatureFromDensityInternalEnergy( + const Real T = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k, j, i), vprim(b, gas::prim::sie(n), k, j, i)); - const Real Tm = eos_d.TemperatureFromDensityInternalEnergy( + const Real Tm = eos_d(n).TemperatureFromDensityInternalEnergy( vprim(b, gas::prim::density(n), k - 1, j, i), vprim(b, gas::prim::sie(n), k - 1, j, i)); diff --git a/src/utils/eos/eos.hpp b/src/utils/eos/eos.hpp index cd81ca7f..c1be3432 100644 --- a/src/utils/eos/eos.hpp +++ b/src/utils/eos/eos.hpp @@ -23,22 +23,15 @@ namespace ArtemisUtils { -// Maximum size of lambda array for optional extra EOS arguments. -// As it happens this must be >= 1 for device code. -static constexpr int lambda_max_vals = 1; - // Variant containing all EOSs to be used in Artemis. using EOS = singularity::Variant #ifdef SPINER_USE_HDF , - singularity::UnitSystem -#ifdef WITH_SESAME - , - singularity::UnitSystem, - singularity::UnitSystem -#endif + singularity::UnitSystem, + singularity::UnitSystem, + singularity::UnitSystem #endif >; diff --git a/src/utils/fluxes/fluid_fluxes.hpp b/src/utils/fluxes/fluid_fluxes.hpp index 8d35c205..e003bbce 100644 --- a/src/utils/fluxes/fluid_fluxes.hpp +++ b/src/utils/fluxes/fluid_fluxes.hpp @@ -87,9 +87,9 @@ TaskStatus CalculateFluxesImpl(MeshData *md, PKG &pkg, PRIM vp, FLUX vflx, const bool three_d = (pm->ndim > 2); // Adiabatic index, if used - EOS eos; + ParArray1D eos; if constexpr (F == Fluid::gas) { - eos = pkg->template Param("eos_d"); + eos = pkg->template Param>("eos_d"); } const auto &cpars = diff --git a/src/utils/fluxes/reconstruction/reconstruction.hpp b/src/utils/fluxes/reconstruction/reconstruction.hpp index c1a7011b..4c842e73 100644 --- a/src/utils/fluxes/reconstruction/reconstruction.hpp +++ b/src/utils/fluxes/reconstruction/reconstruction.hpp @@ -49,7 +49,7 @@ struct ReconGradient { //! fluxes when near ~round-off template KOKKOS_INLINE_FUNCTION void -post_recon(const EOS &eos, const Real dfloor, const Real siefloor, +post_recon(const ParArray1D &eos, const Real dfloor, const Real siefloor, parthenon::team_mbr_t const &member, const int dir, const int b, const int k, const int j, const int il, const int iu, const V &q, parthenon::ScratchPad2D &ql, parthenon::ScratchPad2D &qr) { @@ -78,34 +78,34 @@ post_recon(const EOS &eos, const Real dfloor, const Real siefloor, const int IPR = nspecies * 4 + n; const int ISE = nspecies * 5 + n; const int IBL = nspecies * 6 + n; - parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, il, iu, - [&](const int i) { - const int ipl = i + (dir == 1); - Real &dL = ql(IDN, ipl); - Real &pL = ql(IPR, ipl); - Real &eL = ql(ISE, ipl); - Real &bL = qr(IBL, ipl); - Real &dR = qr(IDN, i); - Real &pR = qr(IPR, i); - Real &eR = qr(ISE, i); - Real &bR = qr(IBL, i); + parthenon::par_for_inner( + DEFAULT_INNER_LOOP_PATTERN, member, il, iu, [&](const int i) { + const int ipl = i + (dir == 1); + Real &dL = ql(IDN, ipl); + Real &pL = ql(IPR, ipl); + Real &eL = ql(ISE, ipl); + Real &bL = qr(IBL, ipl); + Real &dR = qr(IDN, i); + Real &pR = qr(IPR, i); + Real &eR = qr(ISE, i); + Real &bR = qr(IBL, i); - // Floor everything - dL = std::max(dL, dfloor); - dR = std::max(dR, dfloor); - eL = std::max(eL, siefloor); - eR = std::max(eR, siefloor); + // Floor everything + dL = std::max(dL, dfloor); + dR = std::max(dR, dfloor); + eL = std::max(eL, siefloor); + eR = std::max(eR, siefloor); - // Only correct these if something is wrong - if ((pL <= 0.0) || (bL <= 0.0)) { - pL = eos.PressureFromDensityInternalEnergy(dL, eL); - bL = eos.BulkModulusFromDensityInternalEnergy(dL, eL); - } - if ((pR <= 0.0) || (bR <= 0.0)) { - pR = eos.PressureFromDensityInternalEnergy(dR, eR); - bR = eos.BulkModulusFromDensityInternalEnergy(dR, eR); - } - }); + // Only correct these if something is wrong + if ((pL <= 0.0) || (bL <= 0.0)) { + pL = eos(n).PressureFromDensityInternalEnergy(dL, eL); + bL = eos(n).BulkModulusFromDensityInternalEnergy(dL, eL); + } + if ((pR <= 0.0) || (bR <= 0.0)) { + pR = eos(n).PressureFromDensityInternalEnergy(dR, eR); + bR = eos(n).BulkModulusFromDensityInternalEnergy(dR, eR); + } + }); } } } diff --git a/src/utils/fluxes/riemann/hllc.hpp b/src/utils/fluxes/riemann/hllc.hpp index 9cd922f9..619ffced 100644 --- a/src/utils/fluxes/riemann/hllc.hpp +++ b/src/utils/fluxes/riemann/hllc.hpp @@ -47,13 +47,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); @@ -92,8 +92,10 @@ struct RiemannSolver struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/hlle.hpp b/src/utils/fluxes/riemann/hlle.hpp index 0e554747..ae568f93 100644 --- a/src/utils/fluxes/riemann/hlle.hpp +++ b/src/utils/fluxes/riemann/hlle.hpp @@ -53,13 +53,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); @@ -230,13 +230,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/llf.hpp b/src/utils/fluxes/riemann/llf.hpp index 0ae2a0fd..b467355e 100644 --- a/src/utils/fluxes/riemann/llf.hpp +++ b/src/utils/fluxes/riemann/llf.hpp @@ -44,13 +44,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction @@ -175,13 +175,13 @@ template struct RiemannSolver> { template - KOKKOS_INLINE_FUNCTION void operator()(const EOS &eos, const Real c, const Real chat, - parthenon::team_mbr_t const &member, const int b, - const int k, const int j, const int il, - const int iu, const int dir, - const parthenon::ScratchPad2D &wl, - const parthenon::ScratchPad2D &wr, - const V1 &p, const V2 &q, const V3 &vf) const { + KOKKOS_INLINE_FUNCTION void + operator()(const ParArray1D &eos, const Real c, const Real chat, + parthenon::team_mbr_t const &member, const int b, const int k, const int j, + const int il, const int iu, const int dir, + const parthenon::ScratchPad2D &wl, + const parthenon::ScratchPad2D &wr, const V1 &p, const V2 &q, + const V3 &vf) const { using TE = parthenon::TopologicalElement; // Check sensibility of flux direction PARTHENON_REQUIRE(dir > 0 && dir <= 3, "Invalid flux direction!"); diff --git a/src/utils/fluxes/riemann/riemann.hpp b/src/utils/fluxes/riemann/riemann.hpp index 80c42bde..c8ed4022 100644 --- a/src/utils/fluxes/riemann/riemann.hpp +++ b/src/utils/fluxes/riemann/riemann.hpp @@ -26,7 +26,7 @@ template struct RiemannSolver { template KOKKOS_INLINE_FUNCTION void - operator()(const EOS &eos, const Real c, const Real chat, + operator()(const ParArray1D &eos, const Real c, const Real chat, parthenon::team_mbr_t const &member, const int b, const int k, const int j, const int il, const int iu, const int dir, const parthenon::ScratchPad2D &wl, diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt index 2704e440..20244afc 100644 --- a/tst/unit/CMakeLists.txt +++ b/tst/unit/CMakeLists.txt @@ -45,3 +45,4 @@ include(Catch) add_unit_test(test_ideal_gas test_ideal_gas.cpp) add_unit_test(test_ideal_hhe test_ideal_hhe.cpp) add_unit_test(test_geometry test_geometry.cpp) +add_unit_test(test_collision_integrals test_collision_integrals.cpp) diff --git a/tst/unit/README.md b/tst/unit/README.md new file mode 100644 index 00000000..1a2b388f --- /dev/null +++ b/tst/unit/README.md @@ -0,0 +1,121 @@ +# Artemis Unit Tests + + + +This directory contains unit tests for Artemis using the Catch2 testing framework. + +## Building and Running Tests + +To build with unit tests enabled: + +```bash +cmake --preset cpu-debug # or your preferred preset +cmake --build --preset cpu-debug +``` + +To run all unit tests: + +```bash +cd build +ctest +``` + +Or run tests with verbose output: + +```bash +ctest --verbose +``` + +To run only unit tests (excluding regression tests): + +```bash +ctest -L unit +``` + +To run the regression test suite: + +```bash +ctest -L regression +``` + +To exclude regression tests: + +```bash +ctest -LE regression +``` + +To run a specific unit test directly: + +```bash +./tst/unit/test_ideal_gas +``` + +## Regression Tests + +Artemis includes regression tests that run full simulations. These are automatically integrated into CTest: + +- **CPU builds**: Runs `regression.suite` (serial + parallel tests) +- **GPU builds** (CUDA/HIP): Runs `gpu.suite` (serial + GPU-specific tests) + +The regression tests use `tst/run_tests.py` and can take significantly longer than unit tests. Test suites are defined in `tst/suites/`. + +## Adding New Tests + +1. Create a new `.cpp` file in `tst/unit/` +2. Include Catch2 headers: `#include ` +3. Write your tests using Catch2's `TEST_CASE` macro +4. Add your test to `tst/unit/CMakeLists.txt` using `add_unit_test()` + +Example: +```cpp +#include + +TEST_CASE("My test description", "[tag]") { + REQUIRE(1 + 1 == 2); +} +``` + +## Current Tests + +### test_ideal_gas.cpp +Tests for the ideal gas equation of state (IdealGas EOS): +- Gamma parameter recovery +- Pressure, density, and temperature relationships +- Specific internal energy calculations +- UnitSystem wrapper functionality with different unit systems + +### test_ideal_hhe.cpp +Tests for the hydrogen-helium mixture equation of state (IdealHHe EOS): +- Temperature recovery from density and internal energy +- Pressure calculations +- Internal energy monotonicity with temperature +- Different H-He compositions (solar, H-rich, He-rich) +- UnitSystem wrapper with CGS and scaled units + +### test_geometry.cpp +Tests for coordinate system implementations in `src/geometry/`: +- **Cartesian** (x, y, z): Volume, centroids, face areas, unit scale factors, no metric dependence +- **Spherical 3D** (r, θ, φ): Volume, centroids, face areas +- **Spherical 2D** (r, θ): Axisymmetric volume and centroid calculations +- **Spherical 1D** (r): Spherically symmetric volume, radial centroid, face areas, thin shell approximation +- **Cylindrical** (R, φ, z): Volume, radial centroid, face areas, scale factors +- **Axisymmetric** (R, z, φ): Volume, radial centroid, face areas, scale factors +- **Edge cases**: Poles, axis, thin shells + +Each test validates: +- Volume calculations using analytical integration formulas +- Volume-weighted centroids (x1v, x2v, x3v) +- Face areas (AreaX1, AreaX2, AreaX3) +- Metric scale factors (hx1v, hx2v, hx3v) +- Coordinate dependencies (x1dep, x2dep, x3dep) + +## Test Organization + +Tests should be organized by component: +- `test_ideal_gas.cpp` - Tests for ideal gas equation of state +- Add more tests as needed for other components + +## Catch2 Resources + +- [Catch2 Tutorial](https://github.com/catchorg/Catch2/blob/devel/docs/tutorial.md) +- [Assertion Macros](https://github.com/catchorg/Catch2/blob/devel/docs/assertions.md) diff --git a/tst/unit/test_collision_integrals.cpp b/tst/unit/test_collision_integrals.cpp new file mode 100644 index 00000000..a9f8c701 --- /dev/null +++ b/tst/unit/test_collision_integrals.cpp @@ -0,0 +1,273 @@ +//======================================================================================== +// (C) (or copyright) 2026. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +// This file was created in part or in whole by generative AI + +//! \file test_collision_integrals.cpp +//! \brief Catch2 unit tests for the Chapman-Cowling gas-gas collision backends. + +#include + +#include +#include + +#include "drag/collision_integrals.hpp" + +using Catch::Approx; +using Drag::GasDragModel; +using Drag::CollisionIntegrals::DragCoeff; +using Drag::CollisionIntegrals::ThermalRelaxRate; + +namespace { + +double HardSphereCoeff(const double kb_code, const double m0, const double m1, + const double s0, const double s1, const double rho0, + const double rho1, const double T) { + const double sij = 0.5 * (s0 + s1); + const double mij = m0 * m1 / (m0 + m1); + const double nij0 = rho0 / m0; + const double nij1 = rho1 / m1; + const double vrel = std::sqrt(8.0 * kb_code * T / (M_PI * mij)); + return (4.0 / 3.0) * nij0 * nij1 * mij * M_PI * sij * sij * vrel; +} + +double KineticEnergy(const double rho, const double p) { return 0.5 * p * p / rho; } + +void ApplyTwoFluidImpulse(const double K, const double dt, const double rho0, + const double rho1, double &p0, double &p1) { + const double v0 = p0 / rho0; + const double v1 = p1 / rho1; + const double J = K * dt * (v1 - v0) / (1.0 + K * dt * (1.0 / rho0 + 1.0 / rho1)); + p0 += J; + p1 -= J; +} + +double ApplyThermalExchange(const double H, const double dt, const double C0, + const double C1, const double T0, const double T1) { + return H * dt * (T1 - T0) / (1.0 + H * dt * (1.0 / C0 + 1.0 / C1)); +} + +} // namespace + +TEST_CASE("DragCoeff hard_sphere matches Chapman-Cowling form", + "[collision_integrals][hard_sphere]") { + const double kb_code = 1.0; + const double m0 = 1.0; + const double m1 = 16.0; + const double s0 = 1.0e-3; + const double s1 = 2.0e-3; + const double rho0 = 1.0e-4; + const double rho1 = 3.0e-4; + const double T = 2.0e4; + + const double K = DragCoeff(GasDragModel::hard_sphere, kb_code, m0, m1, s0, s1, 0.0, 0.0, + rho0, rho1, T); + const double K_expected = HardSphereCoeff(kb_code, m0, m1, s0, s1, rho0, rho1, T); + + REQUIRE(K == Approx(K_expected).epsilon(1.0e-12)); +} + +TEST_CASE("DragCoeff remains symmetric and positive", + "[collision_integrals][hard_sphere]") { + const double kb_code = 1.0; + const double m0 = 1.0; + const double m1 = 4.0; + const double s0 = 1.0e-3; + const double s1 = 3.0e-3; + const double rho0 = 1.0e-4; + const double rho1 = 2.0e-4; + const double T = 1.0e4; + + const double K01 = DragCoeff(GasDragModel::hard_sphere, kb_code, m0, m1, s0, s1, 0.0, + 0.0, rho0, rho1, T); + const double K10 = DragCoeff(GasDragModel::hard_sphere, kb_code, m1, m0, s1, s0, 0.0, + 0.0, rho1, rho0, T); + + REQUIRE(K01 >= 0.0); + REQUIRE(K01 == Approx(K10).epsilon(1.0e-12)); +} + +TEST_CASE("DragCoeff hard_sphere keeps expected density sigma and temperature scaling", + "[collision_integrals][hard_sphere]") { + const double kb_code = 1.0; + const double m = 2.0; + const double s = 1.0e-3; + const double rho = 1.0e-4; + const double T = 1.0e4; + + const double K_base = + DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, s, s, 0.0, 0.0, rho, rho, T); + const double K_rho = DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, s, s, 0.0, 0.0, + 2.0 * rho, 2.0 * rho, T); + const double K_sigma = DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, 2.0 * s, + 2.0 * s, 0.0, 0.0, rho, rho, T); + const double K_temp = DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, s, s, 0.0, + 0.0, rho, rho, 4.0 * T); + + REQUIRE(K_rho == Approx(4.0 * K_base).epsilon(1.0e-12)); + REQUIRE(K_sigma == Approx(4.0 * K_base).epsilon(1.0e-12)); + REQUIRE(K_temp == Approx(2.0 * K_base).epsilon(1.0e-12)); +} + +TEST_CASE("Two-fluid implicit drag conserves momentum for unequal densities", + "[collision_integrals][momentum]") { + const double K = 5.0; + const double dt = 0.4; + const double rho0 = 1.0; + const double rho1 = 3.0; + double p0 = rho0 * 7.0; + double p1 = rho1 * -2.0; + + const double p_sum_old = p0 + p1; + const double dv_old = std::abs(p1 / rho1 - p0 / rho0); + ApplyTwoFluidImpulse(K, dt, rho0, rho1, p0, p1); + + REQUIRE(p0 + p1 == Approx(p_sum_old).epsilon(1.0e-12)); + REQUIRE(std::abs(p1 / rho1 - p0 / rho0) < dv_old); +} + +TEST_CASE("Two-fluid implicit drag reaches center-of-mass velocity in strong coupling", + "[collision_integrals][momentum]") { + const double K = 1.0e12; + const double dt = 1.0; + const double rho0 = 1.0; + const double rho1 = 4.0; + double p0 = rho0 * 9.0; + double p1 = rho1 * -1.0; + + ApplyTwoFluidImpulse(K, dt, rho0, rho1, p0, p1); + + const double vcom = (rho0 * 9.0 + rho1 * -1.0) / (rho0 + rho1); + REQUIRE(p0 / rho0 == Approx(vcom).epsilon(1.0e-9)); + REQUIRE(p1 / rho1 == Approx(vcom).epsilon(1.0e-9)); +} + +TEST_CASE("Drag kinetic energy loss becomes positive drift heat", + "[collision_integrals][energy]") { + const double K = 10.0; + const double dt = 0.5; + const double rho0 = 1.0; + const double rho1 = 2.0; + const double C0 = 3.0; + const double C1 = 5.0; + double p0 = rho0 * 8.0; + double p1 = rho1 * -1.0; + double e0 = 20.0; + double e1 = 30.0; + double u0 = 12.0; + double u1 = 18.0; + + const double ke0_old = KineticEnergy(rho0, p0); + const double ke1_old = KineticEnergy(rho1, p1); + const double total_old = e0 + e1; + + ApplyTwoFluidImpulse(K, dt, rho0, rho1, p0, p1); + + const double ke0_new = KineticEnergy(rho0, p0); + const double ke1_new = KineticEnergy(rho1, p1); + const double dke0 = ke0_new - ke0_old; + const double dke1 = ke1_new - ke1_old; + const double qdrag = -(dke0 + dke1); + + REQUIRE(qdrag > 0.0); + + e0 += dke0; + e1 += dke1; + const double q0 = qdrag * C0 / (C0 + C1); + const double q1 = qdrag * C1 / (C0 + C1); + e0 += q0; + e1 += q1; + u0 += q0; + u1 += q1; + + REQUIRE(e0 + e1 == Approx(total_old).epsilon(1.0e-12)); + REQUIRE((u0 - 12.0) + (u1 - 18.0) == Approx(qdrag).epsilon(1.0e-12)); +} + +TEST_CASE("ThermalRelaxRate is positive and symmetric as a pair coefficient", + "[collision_integrals][thermal]") { + const double kb_code = 1.0; + const double m0 = 1.0; + const double m1 = 16.0; + const double s0 = 1.0e-3; + const double s1 = 2.0e-3; + const double e0 = 0.0; + const double e1 = 0.0; + const double dof0 = 3.0; + const double dof1 = 3.0; + const double cv0 = 1.5; + const double cv1 = 1.5 / 16.0; + const double rho0 = 1.0e-4; + const double rho1 = 2.0e-4; + const double T = 1.0e4; + + const double H01 = ThermalRelaxRate(GasDragModel::hard_sphere, kb_code, m0, m1, s0, s1, + e0, e1, dof0, dof1, cv0, cv1, rho0, rho1, T); + const double H10 = ThermalRelaxRate(GasDragModel::hard_sphere, kb_code, m1, m0, s1, s0, + e1, e0, dof1, dof0, cv1, cv0, rho1, rho0, T); + + REQUIRE(H01 >= 0.0); + REQUIRE(H01 == Approx(H10).epsilon(1.0e-12)); +} + +TEST_CASE("Implicit thermal exchange conserves energy and reaches weighted equilibrium", + "[collision_integrals][thermal]") { + const double H = 1.0e12; + const double dt = 1.0; + const double C0 = 2.0; + const double C1 = 5.0; + const double T0 = 12.0; + const double T1 = 4.0; + + const double q = ApplyThermalExchange(H, dt, C0, C1, T0, T1); + const double T0_new = T0 + q / C0; + const double T1_new = T1 - q / C1; + const double Teq = (C0 * T0 + C1 * T1) / (C0 + C1); + + REQUIRE(C0 * T0_new + C1 * T1_new == Approx(C0 * T0 + C1 * T1).epsilon(1.0e-12)); + REQUIRE(T0_new == Approx(Teq).epsilon(1.0e-9)); + REQUIRE(T1_new == Approx(Teq).epsilon(1.0e-9)); +} + +TEST_CASE( + "LJ collision integrals enhance cold coupling and approach hard-sphere at high T*", + "[collision_integrals][lj]") { + const double kb_code = 1.0; + const double m = 1.0; + const double s = 1.0e-3; + const double eps = 1.0e4; + const double dof = 3.0; + const double cv = 1.5; + const double rho = 1.0e-4; + const double T_low = 5.0e3; + const double T_high = 1.0e6; + + const double K_lj_low = + DragCoeff(GasDragModel::lj, kb_code, m, m, s, s, eps, eps, rho, rho, T_low); + const double K_hs_low = DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, s, s, 0.0, + 0.0, rho, rho, T_low); + const double K_lj_high = + DragCoeff(GasDragModel::lj, kb_code, m, m, s, s, eps, eps, rho, rho, T_high); + const double K_hs_high = DragCoeff(GasDragModel::hard_sphere, kb_code, m, m, s, s, 0.0, + 0.0, rho, rho, T_high); + + const double H_lj_low = ThermalRelaxRate(GasDragModel::lj, kb_code, m, m, s, s, eps, + eps, dof, dof, cv, cv, rho, rho, T_low); + const double H_lj_high = ThermalRelaxRate(GasDragModel::lj, kb_code, m, m, s, s, eps, + eps, dof, dof, cv, cv, rho, rho, T_high); + + REQUIRE(K_lj_low / K_hs_low > 1.0); + REQUIRE(K_lj_low / K_hs_low > K_lj_high / K_hs_high); + REQUIRE(H_lj_low / K_lj_low > 0.0); + REQUIRE(H_lj_high / K_lj_high > 0.0); +}