Dust coagulation

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 = dust_coagulation           # name of the pgen
+coordinates = cartesian  # coordinate 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	# tab data dump
+dt        = 62.8 	# time increment between outputs
+
+<parthenon/output2>
+file_type = rst
+dt 	  = 628
+
+<parthenon/time>
+nlim       = -1   	# cycle limit
+tlim       = 6280. 	# time limit
+integrator = rk2    	# time integration algorithm
+ncycle_out = 1    	# interval for stdout summary info
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+# numlevel = 1
+
+nx1    = 16  		# Number of zones in X1-direction
+x1min  = 1.35e15   	# minimum value of X1
+x1max  = 1.64e15  	# 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
+gamma = 1.00001
+iso_sound_speed = 0.05	# code unit
+
+<dust>
+cfl = 0.9
+nspecies = 121
+size_input = logspace
+
+scr_level = 1		# for reconstruction
+
+surface_density_flag = true  # surface or volume density 
+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_adaptiveStep = false
+coag_mom_preserve     = false
+coag_info_out	      = true
+
+<problem>
+nInit_dust = 11		# number of dust species initially
+dust_to_gas = 0.01
+nstep1Coag = 10
+mstar = 1.0		# central star masss (m_sun)
+r0_length = 10.0	# 1 code unit (AU)
+rho0 = 2e-5		# code unit value of 1 m_star/r0_length^2
+const_coag_omega = true # using omega at r0_length
+

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 = dust_coagulation           # name of the pgen
+coordinates = cartesian  # coordinate 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	# tab data dump
+dt        = 62.8 	# time increment between outputs
+
+<parthenon/output2>
+file_type = rst
+dt 	  = 628
+
+<parthenon/time>
+nlim       = 5   	# cycle limit
+tlim       = 6280. 	# time limit
+integrator = rk2    	# time integration algorithm
+ncycle_out = 1    	# interval for stdout summary info
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+# numlevel = 1
+
+nx1    = 16  		# Number of zones in X1-direction
+x1min  = 1.35e15   	# minimum value of X1
+x1max  = 1.64e15  	# 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
+gamma = 1.00001
+iso_sound_speed = 0.05	# code unit
+
+<dust>
+cfl = 0.9
+nspecies = 121
+size_input = logspace
+
+scr_level = 1		# for reconstruction
+
+surface_density_flag = false  # surface or volume density 
+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_adaptiveStep = false
+coag_mom_preserve     = false
+coag_info_out	      = true
+
+<problem>
+nInit_dust = 11		# number of dust species initially
+dust_to_gas = 0.01
+nstep1Coag = 1
+mstar = 1.0		# central star masss (m_sun)
+r0_length = 10.0	# 1 code unit (AU)
+rho0 = 2e-5		# code unit value of 1 m_star/r0_length^2
+const_coag_omega = true # using omega at r0_length
+

inputs/dust/dust_collision.in

@@ -0,0 +1,91 @@
+# ========================================================================================
+#  (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 = dust_collision           # name of the pgen
+coordinates = cartesian  # coordinate system
+
+<parthenon/job>
+problem_id = dc_B  # problem ID: basename of output filenames
+
+<parthenon/output1>
+variables = gas.prim.density,   &
+            gas.prim.velocity,  &
+            gas.prim.pressure,  &
+            dust.prim.density,  &
+            dust.prim.velocity
+file_type = hdf5     # tab data dump
+dt        = 0.05  # time increment between outputs
+
+<parthenon/time>
+nlim       = -1     # cycle limit
+tlim       = 0.05   # time limit
+integrator = rk2    # time integration algorithm
+ncycle_out = 1      # interval for stdout summary info
+dt_user    = 0.005
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+# numlevel = 1
+
+nx1    = 16        # Number of zones in X1-direction
+x1min  = 0.0       # minimum value of X1
+x1max  = 1.0       # maximum value of X1
+ix1_bc = periodic  # Inner-X1 boundary condition flag
+ox1_bc = periodic  # Outer-X1 boundary condition flag
+
+nx2    = 1         # Number of zones in X2-direction
+x2min  = -0.5      # minimum value of X2
+x2max  =  0.5      # 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
+drag = true
+
+<gas>
+cfl = 0.9
+gamma = 1.4
+iso_sound_speed = 0.1
+
+<dust>
+cfl = 0.9
+nspecies = 2
+
+sizes = 1.0, 1.0
+
+<dust/stopping_time>
+type = constant
+tau = 0.01, 0.001  # for iprob = 1 
+#tau = 2.0, 1.0    # for iprob = 2
+
+<drag>
+type = simple_dust
+
+
+<problem>
+iprob = 1	   # 1 or 2

src/CMakeLists.txt

@@ -27,6 +27,8 @@ set (SRC_LIST
 
   dust/dust.cpp
   dust/dust.hpp
+  dust/coagulation/coagulation.cpp
+  dust/coagulation/coagulation.hpp
 
   gas/gas.cpp
   gas/gas.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"
@@ -30,6 +31,8 @@
 
 namespace artemis {
 
+std::vector<TaskCollectionFnPtr> OperatorSplitTasks;
+
 //----------------------------------------------------------------------------------------
 //! \fn  Packages_t Artemis::ProcessPackages
 //! \brief Process and initialize relevant packages
@@ -63,6 +66,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);
   artemis->AddParam("do_gas", do_gas);
   artemis->AddParam("do_dust", do_dust);
   artemis->AddParam("do_gravity", do_gravity);
@@ -74,6 +78,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);
   PARTHENON_REQUIRE(!(do_cooling) || (do_cooling && do_gas),
                     "Cooling requires the gas package, but there is not gas!");
   PARTHENON_REQUIRE(!(do_viscosity) || (do_viscosity && do_gas),
@@ -82,6 +87,8 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
                     "Conduction requires the gas package, but there is not gas!");
   PARTHENON_REQUIRE(!(do_radiation) || (do_radiation && do_gas),
                     "Radiation requires the gas package, but there is not gas!");
+  PARTHENON_REQUIRE(!(do_coagulation) || (do_coagulation && do_dust),
+                    "Coagulation requires the dust package, but there is not dust!");
 
   // Set coordinate system
   const int ndim = ProblemDimension(pin.get());
@@ -98,6 +105,10 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get()));
   if (do_drag) packages.Add(Drag::Initialize(pin.get()));
   if (do_nbody) packages.Add(NBody::Initialize(pin.get()));
+  if (do_coagulation) {
+    auto &dustPars = packages.Get("dust")->AllParams();
+    packages.Add(Dust::Coagulation::Initialize(pin.get(), dustPars));
+  }
   if (do_radiation) {
     auto eos_h = packages.Get("gas")->Param<EOS>("eos_h");
     auto opacity_h = packages.Get("gas")->Param<Opacity>("opacity_h");

src/artemis.hpp

@@ -71,7 +71,7 @@ ARTEMIS_VARIABLE(dust.prim, velocity);
 #undef ARTEMIS_VARIABLE
 
 // TaskCollection function pointer for operator split tasks
-using TaskCollectionFnPtr = TaskCollection (*)(Mesh *pm, const Real time, const Real dt);
+using TaskCollectionFnPtr = TaskCollection (*)(Mesh *pm, parthenon::SimTime &tm);
 
 // Constants that enumerate...
 // ...Coordinate systems
@@ -148,6 +148,7 @@ inline int ProblemDimension(parthenon::ParameterInput *pin) {
 
 namespace artemis {
 extern std::function<AmrTag(MeshBlockData<Real> *mbd)> ProblemCheckRefinementBlock;
+// extern std::vector<TaskCollectionFnPtr> OperatorSplitTasks;
 } // namespace artemis
 
 #endif // ARTEMIS_ARTEMIS_HPP_