Dust Coagulation v2

CMakeLists.txt

@@ -148,6 +148,7 @@ if (NOT TARGET parthenon)
   set(PARTHENON_ENABLE_INIT_PACKING ON CACHE BOOL "" FORCE)
   set(PARTHENON_LINT_DEFAULT OFF CACHE BOOL "" FORCE)
   set(PARTHENON_DISABLE_EXAMPLES ON CACHE BOOL "" FORCE)
+  set(PARTHENON_DISABLE_HDF5_COMPRESSION "Disable HDF5 compression" ON)
   set(BUILD_TESTING OFF CACHE BOOL "" FORCE)
   add_subdirectory(external/parthenon parthenon)
 else()

CMakePresets.json

@@ -64,8 +64,7 @@
         "CMAKE_BUILD_TYPE": "Debug",
         "Kokkos_ENABLE_DEBUG_BOUNDS_CHECK": "ON",
         "CMAKE_CXX_COMPILER": "CC",
-        "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON"
+        "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include"
       }
     },
     {
@@ -74,8 +73,7 @@
         "CMAKE_MAKE_PROGRAM": "$env{MAKE_PROGRAM}",
         "CMAKE_BUILD_TYPE": "RelWithDebInfo",
         "CMAKE_CXX_COMPILER": "CC",
-        "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON"
+        "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include"
       }
     },
     {
@@ -89,7 +87,6 @@
         "ARTEMIS_ENABLE_CUDA": "ON",
         "Kokkos_ARCH_AMPERE80": "ON",
         "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON",
         "NUM_GPU_DEVICES_PER_NODE": "4"
       }
     },
@@ -103,7 +100,6 @@
         "ARTEMIS_ENABLE_CUDA": "ON",
         "Kokkos_ARCH_AMPERE80": "ON",
         "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON",
         "NUM_GPU_DEVICES_PER_NODE": "4"
       }
     },
@@ -118,7 +114,6 @@
         "ARTEMIS_ENABLE_CUDA": "ON",
         "Kokkos_ARCH_HOPPER90": "ON",
         "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON",
         "NUM_GPU_DEVICES_PER_NODE": "4"
       }
     },
@@ -131,7 +126,6 @@
         "CMAKE_CXX_COMPILER": "$env{ARTEMIS_HOME}/external/parthenon/external/Kokkos/bin/nvcc_wrapper",
         "ARTEMIS_ENABLE_CUDA": "ON",
         "HDF5_INCLUDE_DIR": "$env{HDF5_ROOT}/include",
-        "PARTHENON_DISABLE_HDF5_COMPRESSION": "ON",
         "Kokkos_ARCH_HOPPER90": "ON",
         "NUM_GPU_DEVICES_PER_NODE": "4"
       }

inputs/dust/dust_coagulation.in

@@ -0,0 +1,107 @@
+# ========================================================================================
+#  (C) (or copyright) 2023-2024. 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.
+# ========================================================================================
+
+<artemis>
+problem = coag               # name of the pgen
+coordinates = cartesian      # coordinate system
+length = 1.495978707e14      # 10AU in cgs unit
+mass = 1.988416e33           # central-star mass (M_sun) in cgs unit
+time = 1.5871820329271033e8  # cgs unit in year/(2*pi)*(r0**3/mstar)**0.5
+temperature = 2.4753531e4    # kevin (vel**2/kb*mu*mp)
+physical_units = cgs         # physical unit system
+
+<parthenon/job>
+problem_id = coag  # problem ID: basename of output filenames
+
+<parthenon/output0>
+file_type = hst
+dt = 0.628
+
+<parthenon/output1>
+variables = gas.prim.density, &
+            dust.prim.density
+file_type = hdf5  # hdf5 data dump
+dt = 62.8         # time increment between outputs
+
+<parthenon/output2>
+file_type = rst
+dt = 628.0
+
+<parthenon/time>
+nlim = -1         # cycle limit
+tlim = 3768.0     # time limit
+integrator = rk2  # time integration algorithm
+ncycle_out = 10   # interval for stdout summary info
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+
+nx1    = 16        # Number of zones in X1-direction
+x1min  = 9.0       # minimum value of X1
+x1max  = 11.0      # maximum value of X1
+ix1_bc = outflow   # Inner-X1 boundary condition flag
+ox1_bc = outflow   # Outer-X1 boundary condition flag
+
+nx2    = 1         # Number of zones in X2-direction
+x2min  = 0.0       # minimum value of X2
+x2max  = 6.283185307179586  # maximum value of X2
+ix2_bc = periodic  # Inner-X2 boundary condition flag
+ox2_bc = periodic  # Outer-X2 boundary condition flag
+
+nx3    = 1         # Number of zones in X3-direction
+x3min  = -0.5      # minimum value of X3
+x3max  =  0.5      # maximum value of X3
+ix3_bc = periodic  # Inner-X3 boundary condition flag
+ox3_bc = periodic  # Outer-X3 boundary condition flag
+
+<parthenon/meshblock>
+nx1 = 16  # Number of cells in each MeshBlock, X1-dir
+nx2 = 1   # Number of cells in each MeshBlock, X2-dir
+nx3 = 1   # Number of cells in each MeshBlock, X3-dir
+
+<physics>
+gas = true
+dust = true
+coagulation = true
+
+<gas>
+cfl = 0.9        # gas cfl number
+gamma = 1.00001  # adiabatic index
+mu = 2.3         # mean molecular weight
+
+<dust> 
+cfl = 0.9              # dust cfl number
+dfloor = 2.0e-25       # density floor=1e-20*rho0
+nspecies = 121         # number of dust size
+size_input = logspace  # dust size distribution
+scr_level = 1          # for hier. scratch registers
+grain_density = 1.25   # dust internal density g/cc
+min_size = 1e-5        # cm
+max_size = 10.0        # cm
+
+<dust/coagulation>
+vfrag = 1000.0  # cm/s
+coag_int = 1
+coag_use_adaptive_step = false
+coag_mom_preserve = false
+coag_info_out = true
+nstep_coag = 10          # number of steps for one coagulation call
+const_coag_omega = true  # using omega at r0_length
+surface_density_flag = true  # surface or volume density 
+
+<problem>
+h0 = 0.05           # scaled height at r0
+rho0 = 2.0e-5       # rh0 = sigma_0/sqrt(2*pi)/h0
+ninit_dust = 11     # number of dust species initially
+dust_to_gas = 0.01  # init dust-to-gas ratio

inputs/dust/dust_coagulation_den.in

@@ -0,0 +1,107 @@
+# ========================================================================================
+#  (C) (or copyright) 2023-2024. 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.
+# ========================================================================================
+
+<artemis>
+problem = coag               # name of the pgen
+coordinates = cartesian      # coordinate system
+length = 1.495978707e14      # 10AU in cgs unit
+mass = 1.988416e33           # central-star mass (M_sun) in cgs unit
+time = 1.5871820329271033e8  # cgs unit in year/(2*pi)*(r0**3/mstar)**0.5
+temperature = 2.4753531e4    # kevin (vel**2/kb*mu*mp)
+physical_units = cgs         # physical unit system
+
+<parthenon/job>
+problem_id = coag  # problem ID: basename of output filenames
+
+<parthenon/output0>
+file_type = hst
+dt = 0.628
+
+<parthenon/output1>
+variables = gas.prim.density, &
+            dust.prim.density
+file_type = hdf5  # hdf5 data dump
+dt = 62.8         # time increment between outputs
+
+<parthenon/output2>
+file_type = rst
+dt = 628.0
+
+<parthenon/time>
+nlim = -1         # cycle limit
+tlim = 3768.0     # time limit
+integrator = rk2  # time integration algorithm
+ncycle_out = 10   # interval for stdout summary info
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+
+nx1    = 16        # Number of zones in X1-direction
+x1min  = 9.0       # minimum value of X1
+x1max  = 11.0      # maximum value of X1
+ix1_bc = outflow   # Inner-X1 boundary condition flag
+ox1_bc = outflow   # Outer-X1 boundary condition flag
+
+nx2    = 1         # Number of zones in X2-direction
+x2min  = 0.0       # minimum value of X2
+x2max  = 6.283185307179586  # maximum value of X2
+ix2_bc = periodic  # Inner-X2 boundary condition flag
+ox2_bc = periodic  # Outer-X2 boundary condition flag
+
+nx3    = 1         # Number of zones in X3-direction
+x3min  = -0.5      # minimum value of X3
+x3max  =  0.5      # maximum value of X3
+ix3_bc = periodic  # Inner-X3 boundary condition flag
+ox3_bc = periodic  # Outer-X3 boundary condition flag
+
+<parthenon/meshblock>
+nx1 = 16  # Number of cells in each MeshBlock, X1-dir
+nx2 = 1   # Number of cells in each MeshBlock, X2-dir
+nx3 = 1   # Number of cells in each MeshBlock, X3-dir
+
+<physics>
+gas = true
+dust = true
+coagulation = true
+
+<gas>
+cfl = 0.9        # gas cfl number
+gamma = 1.00001  # adiabatic index
+mu = 2.3         # mean molecular weight
+
+<dust> 
+cfl = 0.9              # dust cfl number
+dfloor = 1.59577e-24   # density floor=1e-20*rho0
+nspecies = 121         # number of dust size
+size_input = logspace  # dust size distribution
+scr_level = 1          # for hier. scratch registers
+grain_density = 1.25   # dust internal density g/cc
+min_size = 1e-5        # cm
+max_size = 10.0        # cm
+
+<dust/coagulation>
+vfrag = 1000.0  # cm/s
+coag_int = 1
+coag_use_adaptive_step = false
+coag_mom_preserve = false
+coag_info_out = true
+nstep_coag = 10          # number of steps for one coagulation call
+const_coag_omega = true  # using omega at r0_length
+surface_density_flag = false  # surface or volume density 
+
+<problem>
+h0 = 0.05           # scaled height at r0
+rho0 = 1.59577e-4   # rh0 = sigma_0/sqrt(2*pi)/h0
+ninit_dust = 11     # number of dust species initially
+dust_to_gas = 0.01  # init dust-to-gas ratio

src/CMakeLists.txt

@@ -26,6 +26,8 @@ set (SRC_LIST
 
   dust/dust.cpp
   dust/dust.hpp
+  dust/coagulation/coagulation.cpp
+  dust/coagulation/coagulation.hpp
 
   gas/gas.cpp
   gas/gas.hpp
@@ -57,6 +59,7 @@ set (SRC_LIST
   pgen/advection.hpp
   pgen/beam.hpp
   pgen/blast.hpp
+  pgen/coag.hpp
   pgen/conduction.hpp
   pgen/constant.hpp
   pgen/disk.hpp

src/artemis.cpp

@@ -15,6 +15,7 @@
 #include "artemis.hpp"
 #include "artemis_driver.hpp"
 #include "drag/drag.hpp"
+#include "dust/coagulation/coagulation.hpp"
 #include "dust/dust.hpp"
 #include "gas/cooling/cooling.hpp"
 #include "gas/gas.hpp"
@@ -96,6 +97,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   const bool do_viscosity = pin->GetOrAddBoolean("physics", "viscosity", false);
   const bool do_conduction = pin->GetOrAddBoolean("physics", "conduction", false);
   const bool do_radiation = pin->GetOrAddBoolean("physics", "radiation", false);
+  const bool do_coagulation = pin->GetOrAddBoolean("physics", "coagulation", false);
 
   // Determine input file specified algorithms
   const bool do_imc = do_radiation && pin->DoesBlockExist("radiation/imc");
@@ -114,6 +116,8 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
                     "Radiation requires the gas package, but there is not gas!");
   PARTHENON_REQUIRE(!(do_imc && do_moment),
                     "Cannot simultaneously evolve IMC and moments radiation");
+  PARTHENON_REQUIRE(!(do_coagulation) || (do_coagulation && do_dust),
+                    "Coagulation requires the dust package, but there is not dust!");
 
   // Store configuration choices in params
   artemis->AddParam("do_gas", do_gas);
@@ -127,6 +131,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   artemis->AddParam("do_conduction", do_conduction);
   artemis->AddParam("do_diffusion", do_conduction || do_viscosity);
   artemis->AddParam("do_radiation", do_radiation);
+  artemis->AddParam("do_coagulation", do_coagulation);
   artemis->AddParam("do_imc", do_imc);
   artemis->AddParam("do_moment", do_moment);
   artemis->AddParam("do_shear", do_shear);
@@ -149,6 +154,14 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   if (do_rotating_frame) 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()));
+
+  // Operator split dust coagulation
+  if (do_coagulation) {
+    auto &dust_params = packages.Get("dust")->AllParams();
+    packages.Add(Dust::Coagulation::Initialize(pin.get(), dust_params, units, constants));
+  }
+
+  // Operator split radiation
   if (do_radiation) {
     // Top-level radiation package
     packages.Add(Radiation::Initialize(pin.get(), constants, do_imc));

src/artemis.hpp

@@ -67,6 +67,9 @@ namespace prim {
 ARTEMIS_VARIABLE(dust.prim, density);
 ARTEMIS_VARIABLE(dust.prim, velocity);
 } // namespace prim
+namespace coag {
+ARTEMIS_VARIABLE(dust.coag, ncalls);
+}
 } // namespace dust
 
 namespace rad {

src/artemis_driver.cpp

@@ -22,6 +22,7 @@
 #include "artemis.hpp"
 #include "artemis_driver.hpp"
 #include "drag/drag.hpp"
+#include "dust/coagulation/coagulation.hpp"
 #include "dust/dust.hpp"
 #include "gas/cooling/cooling.hpp"
 #include "gas/gas.hpp"
@@ -69,6 +70,7 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in
   do_diffusion = do_viscosity || do_conduction;
   do_imc = artemis_pkg->template Param<bool>("do_imc");
   do_moment = artemis_pkg->template Param<bool>("do_moment");
+  do_coagulation = artemis_pkg->template Param<bool>("do_coagulation");
 
   // Moments integrator
   if (do_moment) {
@@ -131,6 +133,10 @@ TaskListStatus ArtemisDriver<GEOM>::Step() {
   if (do_moment) status = Moments::MomentsDriver<GEOM>(pmesh, tm, rad_integrator.get());
   if (status != TaskListStatus::complete) return status;
 
+  // Operator split, dust coagulation
+  if (do_coagulation) status = Dust::Coagulation::CoagulationDriver<GEOM>(pmesh, tm);
+  if (status != TaskListStatus::complete) return status;
+
   // Compute new dt, (de)refine, and handle sparse (if enabled)
   status = PostStepTasks().Execute();
 

src/artemis_driver.hpp

@@ -54,7 +54,7 @@ class ArtemisDriver : public EvolutionDriver {
   StateDescriptor *artemis_pkg;
   bool do_gas, do_dust, do_moment, do_imc;
   bool do_gravity, do_nbody, do_rotating_frame, do_shear;
-  bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion;
+  bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion, do_coagulation;
   const bool is_restart;
 };