Dust coagulation #25

shengtai · 2024-12-04T16:16:59Z

Background

I merged develop branch with my dust module (feedback + coagulation)

Description of Changes

change in the dust.cpp, dust.hpp, artemis.hpp artemis.cpp, artemis_driver.cpp,
added new files for coagulation,
added new input files for dust testing

Checklist

New features are documented
Tests added for bug fixes and new features
(@lanl.gov employees) Update copyright on changed files

adamdempsey90 · 2024-12-04T19:11:08Z

you need to remove all of your dust drag code from this branch

brryan · 2024-12-05T23:01:42Z

This PR also needs a CI-enrolled test associated with the dust_coagulation.in input; see e.g. https://github.com/lanl/artemis/blob/develop/tst/scripts/drag/drag.py for an example test script. Let us know on teams if you want help with working the test into the CI so it runs automatically and checks against a known answer

adamdempsey90 · 2024-12-06T19:32:05Z

src/artemis.hpp

@@ -67,11 +67,13 @@ namespace prim {
 ARTEMIS_VARIABLE(dust.prim, density);
 ARTEMIS_VARIABLE(dust.prim, velocity);
 } // namespace prim
+ARTEMIS_VARIABLE(dust, stopping_time);


This needs to be removed

It was removed in the new updated code

adamdempsey90 · 2024-12-06T19:33:22Z

src/artemis_driver.cpp

+  // Execute operator split physics
+  for (auto &fn : OperatorSplitTasks) {
+    status = fn(pmesh, tm, integrator->dt).Execute();
+    if (status != TaskListStatus::complete) return status;
+  }
+


I would prefer that instead of all this OperatorSplitTasks boilerplate throughout the code, we just add the CoagulationTasks directly into the artemis driver. That will be much cleaner and clearer.

Looks like @adamdempsey90 had same thought as my above comment.

I changed as you suggested in the new updated code

adamdempsey90 · 2024-12-06T19:33:53Z

src/artemis_driver.cpp

+        // update dust stopping time and dust diffusivity
+        TaskID dust_stopping_time =
+            tl.AddTask(none, Dust::UpdateDustStoppingTime<GEOM>, u0.get());


This task can be removed

I strongly suggest keeping this task and move all dust stopping time calculation to this routine because this is used in many different places. We can reduce the storage to save only the stopping time only for the first dust species, as other dust species only differ by a factor of dust size.

Stopping time already has a function that calculates it and is used by the drag package, so at the very least any other package that uses stopping time should call that same function. I don't see the benefit of storing this, since calculating it is so cheap. But this is again one of those areas where we first implement something in a potentially non-performant way, evaluate the performance, and then potentially adjust the implementation.

We decided previously to calculate stopping times on the fly, templating on the dust stopping model if necessary, rather than store the extra data. Because the (current) stopping models are simple, it is likely to be more performant to recalculate where necessary rather than retrieve the stored values from memory. Assuming dust stopping is proportional to dust size and hardcoding that into the code is the sort of significant complexity that we would only want to do if we measured later that dust stopping really was a bottleneck.

Still reading through the PR... it sounds like I should expect to see these computed on the fly whenever they are needed... I agree this is what we want.

Regardless, just a warning that by not placing these definitions in FillDerived (which we don't want, probably), there is some ambiguity on what these derived fields hold upon AMR remeshing (e.g., an output of this field may not contain what you expect if the output follows a remeshing event).

I remove this and calculate the stopping time on fly the new code

adamdempsey90 · 2024-12-06T19:36:23Z

src/dust/coagulation/coagulation.cpp

+  const int nm = pin->GetOrAddInteger("dust", "nspecies", 1);
+  cpars.nm = nm;
+  cpars.vfrag = pin->GetOrAddReal("dust", "vfrag", 1.e3); // cm/s
+  cpars.nlim = 1e10;
+  cpars.integrator = pin->GetOrAddInteger("dust", "coag_int", 3);
+  cpars.use_adaptive = pin->GetOrAddBoolean("dust", "coag_use_adaptiveStep", true);
+  cpars.mom_coag = pin->GetOrAddBoolean("dust", "coag_mom_preserve", true);
+  cpars.nCall_mx = pin->GetOrAddInteger("dust", "coag_nsteps_mx", 1000);
+  // dust particle internal density g/cc
+  const Real rho_p = pin->GetOrAddReal("dust", "grain_density", 1.25);


We need to be careful about adding new default values to parameters that have already been read by the dust package. nspecies, grain_density are two examples.

I will modify them to use the same value.

fixed in the new updated code

adamdempsey90 · 2024-12-06T19:37:10Z

src/dust/coagulation/coagulation.cpp

+  const Real s_min = pin->GetReal("dust", "min_size");
+  const Real s_max = pin->GetReal("dust", "max_size");
+  const Real mmin = 4.0 * M_PI / 3.0 * rho_p * std::pow(s_min, 3);
+  const Real mmax = 4.0 * M_PI / 3.0 * rho_p * std::pow(s_max, 3);
+  const Real cond = 1.0 / (1.0 - nm) * std::log(mmin / mmax);
+  const Real conc = std::log(mmin);
+  if (std::exp(cond) > std::sqrt(2.0)) {
+    std::stringstream msg;
+    msg << "### FATAL ERROR in dust with coagulation: using nspecies >"
+        << std::log(mmax / mmin) / (std::log(std::sqrt(2.0))) + 1. << " instead of " << nm
+        << std::endl;
+    PARTHENON_FAIL(msg);
+  }
+
+  ParArray1D<Real> dust_size("dustSize", nm);
+  auto dust_size_host = dust_size.GetHostMirror();
+
+  for (int n = 0; n < nm; ++n) {
+    Real mgrid = std::exp(conc + cond * n);
+    dust_size_host(n) = std::pow(3.0 * mgrid / (4.0 * M_PI * cpars.rho_p), 1. / 3.);
+  }
+  dust_size.DeepCopy(dust_size_host);


The dust package has already read an initialized the dust size array, so this can be removed

This is also related to your question why coagulation must uses the "logspace" in dust size input. I will remove this and only retain the check validation.

fixed in the new updated code

adamdempsey90 · 2024-12-06T19:48:12Z

src/dust/dust.cpp

+    if (parthenon::Globals::my_rank == 0) {
+      for (int n = 0; n < nspecies; n++) {
+        std::cout << "dust_size(" << n << ")=" << h_sizes(n) << std::endl;
+      }
+    }


I don't think we need to print this out

I can remove it

removed it in the new updated code

adamdempsey90 · 2024-12-06T19:48:48Z

src/dust/dust.cpp

+  if (enable_dust_coagulation) {
+    if (cgsunit == NULL) {
+      cgsunit = new CGSUnit;
+    }
+    if (!cgsunit->isSet()) {
+      cgsunit->SetCGSUnit(pin);
+    }
+    params.Add("cgs_unit", cgsunit);
+
+    Real rho_p = rho_p_orig;
+    // density in physical unit correspinding to code unit 1.
+    // const Real rho_g0 = pin->GetOrAddReal("dust", "rho_g0", 1.0);
+    const Real rho_g0 = cgsunit->mass0 / cgsunit->vol0;
+    // re-scale rho_p to account for rho_g0 for stopping_time calculation
+    rho_p /= rho_g0;
+
+    // re-scale rho_p by prefactor coefficient
+    if (cgsunit->isurface_den) {
+      rho_p *= 0.5 * M_PI;
+    } else {
+      rho_p *= std::sqrt(M_PI / 8.);
+      rho_p /= cgsunit->length0;
+    }
+
+    if (parthenon::Globals::my_rank == 0) {
+      std::cout << "rho_p, rho_g0=" << rho_p << " " << rho_g0 << std::endl;
+    }
+    params.Add("rho_p", rho_p);


See earlier comment about not doing unit conversions

fixed in the new committed code

adamdempsey90 · 2024-12-06T19:49:00Z

src/dust/dust.cpp

+  // dust stopping-time
+  m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy,
+                Metadata::Sparse});
+  m.SetSparseThresholds(0.0, 0.0, 0.0);
+  dust->AddSparsePool<dust::stopping_time>(m, control_field, dustids);


Can be removed

see my comment above

removed in the new updated code and calculate the stopping time on fly

adamdempsey90 · 2024-12-06T19:49:10Z

src/dust/dust.cpp

+TaskStatus UpdateDustStoppingTime(MeshData<Real> *md) {
+


This whole function can be removed

removed in the new code

adamdempsey90 · 2024-12-06T19:49:59Z

src/dust/dust.cpp

+  auto userDt = params.template Get<Real>("user_dt");
+  auto nspecies = params.template Get<int>("nspecies");
+
+  const auto cfl_number = params.template Get<Real>("cfl");


There's no reason to modify this function anymore

I can remove it

fixed in the new code

shengtai · 2024-12-06T21:36:59Z

The drag is not false during coagulation test. How then I called your stopping time calculation? From: "Adam M. Dempsey" ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Friday, December 6, 2024 at 2:30 PM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Author ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Sli/dust merge request (PR #25) @adamdempsey90 commented on this pull request.

________________________________ In src/artemis_driver.cpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1874015877__;Iw!!Bt8fGhp8LhKGRg!FHKYjLOh55e8_EtZFGk7hcNMwE9F46DMql6kneCwwDuVGEo3qZ3GN7qfyYLW6ugYP6nV_Tk_DyYOHujb2jZ8kg$>:

+ // update dust stopping time and dust diffusivity

+ TaskID dust_stopping_time = + tl.AddTask(none, Dust::UpdateDustStoppingTime<GEOM>, u0.get()); Stopping time already has a function that calculates it and is used by the drag package, so at the very least any other package that uses stopping time should call that same function. I don't see the benefit of storing this, since calculating it is so cheap. But this is again one of those areas where we first implement something in a potentially non-performant way, evaluate the performance, and then potentially adjust the implementation. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1874015877__;Iw!!Bt8fGhp8LhKGRg!FHKYjLOh55e8_EtZFGk7hcNMwE9F46DMql6kneCwwDuVGEo3qZ3GN7qfyYLW6ugYP6nV_Tk_DyYOHujb2jZ8kg$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZG4QXPWKAV32Q3WWC4D2EIJO5AVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDIOBVHEZTENZZGA__;!!Bt8fGhp8LhKGRg!FHKYjLOh55e8_EtZFGk7hcNMwE9F46DMql6kneCwwDuVGEo3qZ3GN7qfyYLW6ugYP6nV_Tk_DyYOHuhfy3nlVg$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

brryan · 2024-12-09T15:05:24Z

inputs/dust/SI_strat_AB.in

@@ -0,0 +1,135 @@
+# ========================================================================================


Only the two dust coagulation input files show be a part of this PR

this input file is deleted in the new code

brryan · 2024-12-09T15:09:16Z

src/pgen/dust_collision.hpp

+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+
+// NOTE(@pdmullen): The following is taken directly from the open-source


Suggested change

// NOTE(@pdmullen): The following is taken directly from the open-source

// NOTE(@Shengtai): The following is taken directly from the open-source

fixed in the new code

shengtai · 2024-12-09T17:55:33Z

I will remove this in src/artemis_driver.cpp and have a place-holder for Adam’s new routine to calculate stopping time in my coagulation code. From: Ben Ryan ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 9, 2024 at 7:59 AM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Author ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @brryan commented on this pull request.

________________________________ In src/artemis_driver.cpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876129264__;Iw!!Bt8fGhp8LhKGRg!E6TIQ6knxNEWb1UntkXUCHnYaMSwdvzizHDQi_lfV49erpMagaasGz0b_fizL6kSlmX9tdregjfOttoH1fxGEA$>:

+ // update dust stopping time and dust diffusivity

+ TaskID dust_stopping_time = + tl.AddTask(none, Dust::UpdateDustStoppingTime<GEOM>, u0.get()); We decided previously to calculate stopping times on the fly, templating on the dust stopping model if necessary, rather than store the extra data. Because the (current) stopping models are simple, it is likely to be more performant to recalculate where necessary rather than retrieve the stored values from memory. Assuming dust stopping is proportional to dust size and hardcoding that into the code is the sort of significant complexity that we would only want to do if we measured later that dust stopping really was a bottleneck. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876129264__;Iw!!Bt8fGhp8LhKGRg!E6TIQ6knxNEWb1UntkXUCHnYaMSwdvzizHDQi_lfV49erpMagaasGz0b_fizL6kSlmX9tdregjfOttoH1fxGEA$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZG63T4VN7EWNEWA2BST2EWV2TAVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDIOBZGA3DENJTGQ__;!!Bt8fGhp8LhKGRg!E6TIQ6knxNEWb1UntkXUCHnYaMSwdvzizHDQi_lfV49erpMagaasGz0b_fizL6kSlmX9tdregjfOttqXEOTIVw$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

shengtai · 2024-12-09T17:58:01Z

I will move this to coagulation test problem src/pgen/dust_coagulation.hpp. When the artemis unit conversions are ready, I will delete it from the test problem. From: Ben Ryan ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 9, 2024 at 8:00 AM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Author ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @brryan commented on this pull request.

________________________________ In src/dust/coagulation/coagulation.cpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876131360__;Iw!!Bt8fGhp8LhKGRg!BipsjS2KKRz0oVx-CMzvTaRHsak3Fgrjr0mfnuDNLt0arvZTuj0I2a8oGRr34K7WM4xiKPJUYFsYNNJiLcrsFw$>:

+ if (Dust::cgsunit == NULL) {

+ Dust::cgsunit = new CGSUnit; + } + if (!Dust::cgsunit->isSet()) { + Dust::cgsunit->SetCGSUnit(pin); + } Yeah lets ignore any details about units here and just expect a specific unit from the user input file. These dust-specific unit conversions will be replaced by the global #24<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/24__;!!Bt8fGhp8LhKGRg!BipsjS2KKRz0oVx-CMzvTaRHsak3Fgrjr0mfnuDNLt0arvZTuj0I2a8oGRr34K7WM4xiKPJUYFsYNNK_x9LlcQ$>. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876131360__;Iw!!Bt8fGhp8LhKGRg!BipsjS2KKRz0oVx-CMzvTaRHsak3Fgrjr0mfnuDNLt0arvZTuj0I2a8oGRr34K7WM4xiKPJUYFsYNNJiLcrsFw$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZGZ47UFKSZBSXWNVHDD2EWV7FAVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDIOBZGA3DMMBXHE__;!!Bt8fGhp8LhKGRg!BipsjS2KKRz0oVx-CMzvTaRHsak3Fgrjr0mfnuDNLt0arvZTuj0I2a8oGRr34K7WM4xiKPJUYFsYNNLnSAIo5w$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

shengtai · 2024-12-09T17:59:31Z

I have replaced the pow(x,3) with x*x*x. From: Ben Ryan ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 9, 2024 at 8:02 AM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Author ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @brryan commented on this pull request.

________________________________ In src/dust/coagulation/coagulation.cpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876134771__;Iw!!Bt8fGhp8LhKGRg!CRMskvVsk9V6oAdkPDTu95MMzOPlDJ6V9IynKAB0uMN8B0odQJPmI8UHS0lcUh2ZS94VzDAG_VjsbAzaQ8Aa8g$>:

+ ParArray3D<int> cpod_notzero, ParArray3D<Real> cpod_short) {

+ + int ikdelta = coag2DRv::kdelta; + auto kdelta = Kokkos::subview(coag3D, ikdelta, Kokkos::ALL, Kokkos::ALL); + int icoef_fett = coag2DRv::coef_fett; + auto coef_fett = Kokkos::subview(coag3D, icoef_fett, Kokkos::ALL, Kokkos::ALL); + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "initializeCoag1", parthenon::DevExecSpace(), + 0, nm - 1, KOKKOS_LAMBDA(const int i) { + // initialize in_idx(*) array + // initialize kdelta array + for (int j = 0; j < nm; j++) { + kdelta(i, j) = 0.0; + } + kdelta(i, i) = 1.0; + mass_grid(i) = 4.0 * M_PI / 3.0 * rho_p * std::pow(dsize(i), 3); We already have SQR() for std::pow(x,2); if x*x*x shows up in a lot of places we can add a CUBE() { return x*x*x; } function — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876134771__;Iw!!Bt8fGhp8LhKGRg!CRMskvVsk9V6oAdkPDTu95MMzOPlDJ6V9IynKAB0uMN8B0odQJPmI8UHS0lcUh2ZS94VzDAG_VjsbAzaQ8Aa8g$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZGYZY5BRP5U47KEMT3T2EWWHJAVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDIOBZGA3TEMBYGE__;!!Bt8fGhp8LhKGRg!CRMskvVsk9V6oAdkPDTu95MMzOPlDJ6V9IynKAB0uMN8B0odQJPmI8UHS0lcUh2ZS94VzDAG_VjsbAxvO06RVg$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

shengtai · 2024-12-09T18:02:33Z

I will remove “SI_strat_AB.in” as it requires other capacity that does not available in the current artemis and parthenon. I will change “sli/pgen/dust_collision.hpp” as your suggestion. From: Ben Ryan ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 9, 2024 at 8:09 AM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Author ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @brryan commented on this pull request.

________________________________ In inputs/dust/SI_strat_AB.in<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876140718__;Iw!!Bt8fGhp8LhKGRg!D4ZVEl0rqEmoGy-NqFLH0dUbE23CFyqIkeBQtgq2P9NWPmGuxBJJXDND5WvOEwMFDK4e69VPyQ3v0MV6oLHP6A$>:

@@ -0,0 +1,135 @@

+# ======================================================================================== Only the two dust coagulation input files show be a part of this PR

________________________________ In src/pgen/dust_collision.hpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1876150245__;Iw!!Bt8fGhp8LhKGRg!D4ZVEl0rqEmoGy-NqFLH0dUbE23CFyqIkeBQtgq2P9NWPmGuxBJJXDND5WvOEwMFDK4e69VPyQ3v0MUxM-qz1w$>:

+//

+// 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. +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone ***@***.***> and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== + +// ***@***.***): The following is taken directly from the open-source ⬇️ Suggested change -// ***@***.***): The following is taken directly from the open-source +// ***@***.***): The following is taken directly from the open-source — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*pullrequestreview-2489083222__;Iw!!Bt8fGhp8LhKGRg!D4ZVEl0rqEmoGy-NqFLH0dUbE23CFyqIkeBQtgq2P9NWPmGuxBJJXDND5WvOEwMFDK4e69VPyQ3v0MVqTywQ9A$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZG2WG2CVFILINWLBXAD2EWXDLAVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDIOBZGA4DGMRSGI__;!!Bt8fGhp8LhKGRg!D4ZVEl0rqEmoGy-NqFLH0dUbE23CFyqIkeBQtgq2P9NWPmGuxBJJXDND5WvOEwMFDK4e69VPyQ3v0MVoWHrnmA$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

pdmullen

Nice start @shengtai. I think it needs a bit more work to adopt some artemis-isms, particularly as it relates to sparse compatibility and style.

pdmullen · 2024-12-09T19:16:00Z

inputs/dust/SI_strat1_amr_BA.in

+            dust.prim.velocity
+
+file_type = hdf5  # HDF5 data dump
+dt        = 2.0  # time increment between outputs


align comments here and elsewhere.

pdmullen · 2024-12-09T19:16:13Z

inputs/dust/SI_strat1_amr_BA.in

+x2max =  0.02
+level = 1
+
+


remove extraneous whitespace here and elsewhere.

this input is removed in the new code

pdmullen · 2024-12-09T19:19:53Z

src/artemis.cpp

@@ -138,6 +146,25 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
  parthenon::MetadataFlag MetadataOperatorSplit =
      parthenon::Metadata::AddUserFlag("OperatorSplit");

+  if (do_coagulation) {


I feel that coagulation need not use the "anonymous" operator split framework. Even coupling with jaybenne has moved away from this design. I'd mimic the handling of IMC for dust coagulation... e.g.,

if (do_coagulation) status = Dust::OperatorSplitDust<GEOM>(pmesh, tm.time, tm.dt); if (status != TaskListStatus::complete) return status;

maybe we should completely remove the OperatorSplitTasks machinery?

done as you suggested.

pdmullen · 2024-12-09T19:21:05Z

src/artemis.hpp

 } // namespace dust
 #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,


This will break coupling to jaybenne, but I am happy to make the change. If we are passing the whole SimTime object, however, should we still pass dt? Maybe we would want to pass a different dt than tm.dt for something like sub-cycling...

I need SimTime object. It would be better to keep both.

@pdmullen Can subcycling not always be done inside the enrolled task collection? It seems to me to make sense to always provide the SimTime argument here, but maybe I'm missing something.

I would prefer that anything that returns a TaskCollection lives inside the driver class. You then have access to everything in the driver class like SimTime.

pdmullen · 2024-12-09T19:22:29Z

src/artemis_driver.cpp

+  // Execute operator split physics
+  for (auto &fn : OperatorSplitTasks) {
+    status = fn(pmesh, tm, integrator->dt).Execute();
+    if (status != TaskListStatus::complete) return status;
+  }
+


Looks like @adamdempsey90 had same thought as my above comment.

src/dust/dust.cpp

pdmullen · 2024-12-09T19:44:06Z

src/dust/dust.cpp

+  const int nspecies = coag.nm;
+  std::vector<std::string> dust_var_names;
+  for (int n = 0; n < nspecies; n++) {
+    dust_var_names.push_back(dust::prim::density::name() + '_' + std::to_string(n));


Do we really have to do this? I'd like to look at how to clean this up.

If you have simple way to do it, please let me know

Why do you need the names like this? I can't think of a reason why.

adamdempsey90 · 2024-12-10T19:47:52Z

src/dust/dust.cpp

+
+          // for density (g/cc), Stokes number = sqrt(Pi/8)*rho_p*s_p*Omega_k/rho_g/Cs
+          const Real pres = vmesh(b, gas::prim::pressure(0), k, j, i);
+          const Real cs = std::sqrt(gamma * pres / dens_g);


Whenever you need cs, you should get it from the eos via the bulk modulus, e.g., BulkModulusFromDensityInternalEnergy.

done as you suggested

adamdempsey90 · 2024-12-10T19:50:30Z

src/dust/dust.cpp

+          const Real omega1 = Omega_k / time0;
+          const int nm = nspecies;
+
+          ScratchPad1D<Real> rhod(mbr.team_scratch(scr_level), nspecies);


Maybe we should think a bit more about how to store this information. I could see this being okay for when you have > 100 dust species, but maybe not being very performant for < 20 species. Maybe in those smaller cases, you'd want a 2D scratchpad that is (nspecies,nx1) large.

For coagulation run, nspecies is usually larger than 100.

brryan · 2024-12-23T22:12:33Z

src/dust/coagulation/coagulation.cpp

+  CoagParams cpars;
+  const int nm = pin->GetOrAddInteger("dust", "nspecies", 1);
+  cpars.nm = nm;
+  cpars.vfrag = pin->GetOrAddReal("dust", "vfrag", 1.e3); // cm/s


Now that #24 is merged, "unit"-ed quantities like here should be updated to (1) have a specified unit system for input parameters, e.g. cm./s as here (2) converted to code units at the first opportunity, so for example this line should be

cpars.vfrag = pin->GetOrAddReal("dust", "vfrag", 1.e3) * units.GetSpeedPhysicalToCode; // cm/s

where for Initialize functions like this units and constants must be added to the function argument (e.g. if (do_gas) packages.Add(Gas::Initialize(pin.get(), units, constants));), whereas for tasks you can pull the units object out of the artemis package's params.

brryan · 2024-12-23T22:13:57Z

src/dust/coagulation/coagulation.cpp

+  const Real GRAV_CONST = 6.674299999999999e-8; // gravitational const in cm^3 g^-1 s^-2
+  const Real mstar = pin->GetOrAddReal("problem", "mstar", 1.0) * M_SUN;
+  cpars.gm = std::sqrt(GRAV_CONST * mstar);


Suggested change

const Real GRAV_CONST = 6.674299999999999e-8; // gravitational const in cm^3 g^-1 s^-2

const Real mstar = pin->GetOrAddReal("problem", "mstar", 1.0) * M_SUN;

cpars.gm = std::sqrt(GRAV_CONST * mstar);

const Real mstar = pin->GetOrAddReal("problem", "mstar", 1.0) * constants.GetMsolarCode();

cpars.gm = std::sqrt(constants.GetGCode() * mstar);

hardcoded constants need to be replaced with centralized constants from the Constants class introduced in #24

brryan · 2024-12-23T22:14:37Z

src/dust/coagulation/coagulation.cpp

+//! \fn  StateDescriptor Coagulalation::Initialize
+//! \brief Adds intialization function for coagulation package
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Params &dustPars) {
+  auto coag = std::make_shared<StateDescriptor>("coagulation");


Seems to me that this just belongs in the dust package. Thoughts @pdmullen @adamdempsey90?

src/dust/dust.cpp

brryan · 2024-12-23T22:18:02Z

inputs/dust/dust_collision.in

@@ -0,0 +1,91 @@
+# ========================================================================================


This is a pure drag input problem and so is out of scope for this PR unless I'm misunderstanding, so this input file and pgen/dust_collision.hpp should be removed.

shengtai · 2024-12-23T22:54:13Z

When you want to do boundary communication only for part of variables (not all variables), the variable name is needed. Simple “dust::prim:density::name()” does not work for some reason. From: "Adam M. Dempsey" ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 23, 2024 at 3:51 PM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Mention ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @adamdempsey90 commented on this pull request.

________________________________ In src/dust/dust.cpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1896182734__;Iw!!Bt8fGhp8LhKGRg!Avni7ycRIYYET1EQJfh74l8WGXqbLnekRiXrFqDwUKrL0Z3v6VpKvIG6_mgcyaSfF6QhhHAqwhUW4ix3VMVeXA$>:

+

+ const Real time_local = tm.time + dt - (*dtCoag); + const Real dt_local = (*dtCoag); + if (parthenon::Globals::my_rank == 0) { + std::cout << "coagulation at: time,dt,cycle=" << time_local << " " << dt_local << " " + << tm.ncycle << std::endl; + } + + TaskID none(0); + + // Assemble tasks + auto &coag = coag_pkg->template Param<Dust::Coagulation::CoagParams>("coag_pars"); + const int nspecies = coag.nm; + std::vector<std::string> dust_var_names; + for (int n = 0; n < nspecies; n++) { + dust_var_names.push_back(dust::prim::density::name() + '_' + std::to_string(n)); Why do you need the names like this? I can't think of a reason why. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1896182734__;Iw!!Bt8fGhp8LhKGRg!Avni7ycRIYYET1EQJfh74l8WGXqbLnekRiXrFqDwUKrL0Z3v6VpKvIG6_mgcyaSfF6QhhHAqwhUW4ix3VMVeXA$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZG3JMOAG7A6Q3Z6WQ732HCHYRAVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDKMRRGIZTAMZWGM__;!!Bt8fGhp8LhKGRg!Avni7ycRIYYET1EQJfh74l8WGXqbLnekRiXrFqDwUKrL0Z3v6VpKvIG6_mgcyaSfF6QhhHAqwhUW4iycQZkFtA$>. You are receiving this because you were mentioned.Message ID: ***@***.***>

shengtai · 2024-12-23T22:55:09Z

I think this can be removed since both Patrick’s and mine operator split does not use this tastcollection. From: "Adam M. Dempsey" ***@***.***> Reply-To: lanl/artemis ***@***.***> Date: Monday, December 23, 2024 at 3:49 PM To: lanl/artemis ***@***.***> Cc: "Li, Shengtai" ***@***.***>, Mention ***@***.***> Subject: [EXTERNAL] Re: [lanl/artemis] Dust coagulation (PR #25) @adamdempsey90 commented on this pull request.

________________________________ In src/artemis.hpp<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1896181704__;Iw!!Bt8fGhp8LhKGRg!Fwb5yaTmnpgxYodSrLsg3jiHJxAbSbwhUpbo9v1zVfMHz2rqOYuX9pBGTQfBtI9GtRSHBIO4VoL3hlAiL21xUQ$>:

} // namespace dust

#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, I would prefer that anything that returns a TaskCollection lives inside the driver class. You then have access to everything in the driver class like SimTime. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/lanl/artemis/pull/25*discussion_r1896181704__;Iw!!Bt8fGhp8LhKGRg!Fwb5yaTmnpgxYodSrLsg3jiHJxAbSbwhUpbo9v1zVfMHz2rqOYuX9pBGTQfBtI9GtRSHBIO4VoL3hlAiL21xUQ$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AHF3ZG3ETGMX7YIHHZA2LS32HCHQ5AVCNFSM6AAAAABTARH2QOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDKMRRGIZDQOBSGE__;!!Bt8fGhp8LhKGRg!Fwb5yaTmnpgxYodSrLsg3jiHJxAbSbwhUpbo9v1zVfMHz2rqOYuX9pBGTQfBtI9GtRSHBIO4VoL3hlDUxMaPoQ$>. You are receiving this because you were mentioned.Message ID: ***@***.***>

adamdempsey90 · 2024-12-23T23:07:02Z

src/dust/dust.cpp

+TaskListStatus OperatorSplitDust(Mesh *pm, parthenon::SimTime &tm) {
+  auto &dust_pkg = pm->packages.Get("dust");
+  typedef Coordinates C;
+  const C coords = dust_pkg->template Param<Coordinates>("coords");
+
+  if (coords == C::cartesian) {
+    return OperatorSplitDustSelect<C::cartesian>(pm, tm).Execute();
+  } else if (coords == C::spherical1D) {
+    return OperatorSplitDustSelect<C::spherical1D>(pm, tm).Execute();
+  } else if (coords == C::spherical2D) {
+    return OperatorSplitDustSelect<C::spherical2D>(pm, tm).Execute();


OperatorSplitDust is the function that I'm suggesting be put into the driver class. You then have access to SimTime and other things. You can also template it on the coords, so you don't need OperatorSplitDustSelect at all.

shengtai added 2 commits December 3, 2024 23:17

merge with the develop branch

e061712

add files from gitlab/sli/coagulation branch and merge with develop

b191975

shengtai assigned brryan Dec 4, 2024

shengtai added 2 commits December 4, 2024 09:40

format fixed in pgen/dust_coagulation.hpp

63e37b9

Merge branch 'develop' of github.com:lanl/artemis into sli/dust

d82919a

shengtai added 2 commits December 4, 2024 22:21

fixed omega_k calculation in dust.cpp

6ba72b0

remove dust_drag

1beaf24

adamdempsey90 reviewed Dec 6, 2024

View reviewed changes

adamdempsey90 requested a review from pdmullen December 6, 2024 21:39

brryan changed the title ~~Sli/dust merge request~~ Dust coagulation Dec 9, 2024

brryan reviewed Dec 9, 2024

View reviewed changes

pdmullen reviewed Dec 9, 2024

View reviewed changes

adamdempsey90 reviewed Dec 10, 2024

View reviewed changes

shengtai added 5 commits December 11, 2024 13:42

modification based on merge comments

a26caad

Merge branch 'develop' of github.com:lanl/artemis into sli/dust

05da651

Merge branch 'develop' into sli/dust

548a9d1

merge with the head version and modify c2p to p2c in dust.cpp

347ea01

Merge branch 'sli/dust' of github.com:lanl/artemis into sli/dust

11afd99

brryan reviewed Dec 23, 2024

View reviewed changes

adamdempsey90 reviewed Dec 23, 2024

View reviewed changes

		@@ -0,0 +1,135 @@
		# ========================================================================================

	// NOTE(@pdmullen): The following is taken directly from the open-source
	// NOTE(@Shengtai): The following is taken directly from the open-source

		@@ -0,0 +1,91 @@
		# ========================================================================================

		x2max = 0.02
		level = 1

Dust coagulation #25

Are you sure you want to change the base?

Dust coagulation #25

Conversation

shengtai commented Dec 4, 2024

Background

Description of Changes

Checklist

adamdempsey90 commented Dec 4, 2024

brryan commented Dec 5, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shengtai commented Dec 6, 2024 via email

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shengtai commented Dec 9, 2024 via email

shengtai commented Dec 9, 2024 via email

shengtai commented Dec 9, 2024 via email

shengtai commented Dec 9, 2024 via email

pdmullen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shengtai commented Dec 23, 2024 via email

shengtai commented Dec 23, 2024 via email

Choose a reason for hiding this comment