Dust coagulation

inputs/dust/dust_coagulation.in

@@ -0,0 +1,109 @@
+# ========================================================================================
+#  (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
+unit_conversion	= ppd	# ppd (code unit: AU, Msun, year/(2*pi)
+r0_length	= 10.	# AU
+mstar	= 1.0		# M_sun
+rho0	= 2e-5		# prefactor for density profile
+physical_units	= cgs	# 	
+
+<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
+
+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
+iso_sound_speed = 0.05	# code unit
+cv = 1e5		# code unit: 1/(gamma-1)
+
+<dust>
+cfl = 0.9		# dust cfl number
+nspecies = 121		# number of dust size
+size_input = logspace	# dust size distribution
+
+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	# init dust-to-gas ratio
+nstep1Coag = 10		# number of steps for one coagulation call
+const_coag_omega = true # using omega at r0_length
+

inputs/dust/dust_coagulation_den.in

@@ -0,0 +1,109 @@
+# ========================================================================================
+#  (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
+unit_conversion	= ppd	# ppd (code unit: AU, Msun, year/(2*pi)
+r0_length	= 10.	# AU
+mstar	= 1.0		# M_sun
+rho0	= 1.59577e-4	# prefactor for 3D density profile
+physical_units	= cgs	# 	
+
+<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
+
+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
+iso_sound_speed = 0.05	# code unit
+cv = 1e5		# code unit: 1/(gamma-1)
+
+<dust>
+cfl = 0.9		# dust cfl number
+nspecies = 121		# number of dust size
+size_input = logspace	# dust size distribution
+
+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	# init dust-to-gas ratio
+nstep1Coag = 10		# number of steps for one coagulation call
+const_coag_omega = true # using omega at r0_length
+

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"
@@ -31,6 +32,8 @@
 
 namespace artemis {
 
+std::vector<TaskCollectionFnPtr> OperatorSplitTasks;
+
 //----------------------------------------------------------------------------------------
 //! \fn  Packages_t Artemis::ProcessPackages
 //! \brief Process and initialize relevant packages
@@ -70,6 +73,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);
@@ -81,6 +85,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),
@@ -89,6 +94,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());
@@ -105,6 +112,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(), constants));
+  if (do_coagulation) {
+    auto &dustPars = packages.Get("dust")->AllParams();
+    packages.Add(Dust::Coagulation::Initialize(pin.get(), dustPars, units, constants));
+  }
   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_

src/artemis_driver.cpp

@@ -66,6 +66,7 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in
   do_nbody = artemis_pkg->template Param<bool>("do_nbody");
   do_diffusion = do_viscosity || do_conduction;
   do_radiation = artemis_pkg->template Param<bool>("do_radiation");
+  do_coagulation = artemis_pkg->template Param<bool>("do_coagulation");
 
   // NBody initialization tasks
   if (do_nbody) {
@@ -111,6 +112,9 @@ TaskListStatus ArtemisDriver<GEOM>::Step() {
   if (do_radiation) status = IMC::JaybenneIMC<GEOM>(pmesh, tm.time, tm.dt);
   if (status != TaskListStatus::complete) return status;
 
+  if (do_coagulation) status = Dust::OperatorSplitDust<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

@@ -53,11 +53,12 @@ class ArtemisDriver : public EvolutionDriver {
   IntegratorPtr_t integrator, nbody_integrator;
   StateDescriptor *artemis_pkg;
   bool do_gas, do_dust, do_gravity, do_rotating_frame, do_cooling, do_drag, do_viscosity,
-      do_nbody, do_conduction, do_diffusion, do_radiation;
+      do_nbody, do_conduction, do_diffusion, do_radiation, do_coagulation;
   const bool is_restart;
 };
 
-using TaskCollectionFnPtr = TaskCollection (*)(Mesh *pm, const Real time, const Real dt);
+// using TaskCollectionFnPtr = TaskCollection (*)(Mesh *pm, const Real time, const Real
+// dt);
 
 Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin);