Gas-Gas Coupling and Sparse

.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

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

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.
+
+<artemis>
+  problem    = escape_1d
+  coordinates = spherical
+  radial_spacing = logarithmic
+
+<parthenon/job>
+  problem_id = drag
+
+<parthenon/time>
+  nlim      = 1
+  tlim      = 50.0        # ~14 t_dyn (t_dyn = r0/v_esc = 1/sqrt(2))
+  integrator = rk2
+  ncycle_out = 1
+
+<parthenon/mesh>
+  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
+
+<parthenon/meshblock>
+  nx1 = 64
+  nx2 = 1
+  nx3 = 1
+
+<parthenon/output0>
+  file_type  = hdf5
+  dt         = 0.1
+  variables  = gas.prim.density, gas.prim.velocity, gas.prim.sie, gas.prim.temperature, gas.prim.pressure
+
+<parthenon/output1>
+  file_type  = hst
+  dt         = 0.01
+
+# Gas: two-species ideal gas
+<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
+<gravity/point>
+  mass   = 1.0                # G*M in code units
+
+
+<physics>
+  gas = true
+  gravity = true
+  drag = true
+
+# Drag: full Chapman-Cowling coupling
+<drag>
+  type = full
+
+<drag/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
+<problem>
+  # 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)

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.
+
+<artemis>
+  problem    = escape_1d
+  coordinates = spherical
+  radial_spacing = logarithmic
+
+<parthenon/sparse>
+  enable_sparse = true
+  dealloc_count = 1
+  alloc_threshold = 1e-8
+  dealloc_threshold = 1e-10
+
+<parthenon/job>
+  problem_id = drag
+
+<parthenon/time>
+  nlim      = 100
+  tlim      = 50.0        # ~14 t_dyn (t_dyn = r0/v_esc = 1/sqrt(2))
+  integrator = rk2
+  ncycle_out = 1
+
+<parthenon/mesh>
+  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
+
+<parthenon/meshblock>
+  nx1 = 8
+  nx2 = 1
+  nx3 = 1
+
+<parthenon/output0>
+  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
+
+<parthenon/output1>
+  file_type  = hst
+  dt         = 0.01
+
+<gas>
+  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
+<gravity/point>
+  mass   = 1.0                # G*M in code units
+
+
+<physics>
+  gas = true
+  gravity = true
+  drag = true
+
+# Drag: full Chapman-Cowling coupling
+<drag>
+  type = full
+
+<drag/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
+<problem>
+  # 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)

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

src/artemis.cpp

@@ -102,6 +102,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &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<ParameterInput> &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<ParameterInput> &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<ParameterInput> &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>("eos_h");
+      auto eos_h_arr =
+          packages.Get("gas")->Param<parthenon::ParArray1D<ArtemisUtils::EOS>>("eos_h");
       auto opacity_h = packages.Get("gas")->Param<MeanOpacity>("opacity_h");
       auto scattering_h = packages.Get("gas")->Param<MeanScattering>("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<ParameterInput> &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);
 

src/artemis_driver.cpp

@@ -75,11 +75,13 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in
   do_moment = artemis_pkg->template Param<bool>("do_moment");
   do_coagulation = artemis_pkg->template Param<bool>("do_coagulation");
   do_raytrace = artemis_pkg->template Param<bool>("do_raytrace");
+  do_closure = artemis_pkg->template Param<bool>("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<bool>("update_fluxes");
 
+  do_sparse = artemis_pkg->template Param<bool>("do_sparse");
   ndim = pm->ndim;
 
   // Moments integrator
@@ -384,6 +386,10 @@ TaskCollection ArtemisDriver<GEOM>::PostStepTasks() {
     if (pmesh->adaptive) {
       refine = tl.AddTask(new_dt, parthenon::Refinement::Tag<MeshData<Real>>, u0.get());
     }
+    auto dealloc = refine;
+    if (do_sparse) {
+      dealloc = tl.AddTask(refine, parthenon::Update::SparseDealloc, u0.get());
+    }
   }
 
   return tc;

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;
 };
 

src/derived/fill_derived.cpp

@@ -240,12 +240,12 @@ void PrimToCons(T *md) {
   // Extract gas parameters
   Real dflr_gas = Null<Real>();
   Real sieflr_gas = Null<Real>();
-  EOS eos_d;
+  ParArray1D<EOS> eos_d;
   if (do_gas) {
     auto &gas_pkg = pm->packages.Get("gas");
     dflr_gas = gas_pkg->template Param<Real>("dfloor");
     sieflr_gas = gas_pkg->template Param<Real>("siefloor");
-    eos_d = gas_pkg->template Param<EOS>("eos_d");
+    eos_d = gas_pkg->template Param<ParArray1D<EOS>>("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<Real>()};
           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));