Merge branch 'dust_m82' into dust_m82_test

builds/make.type.m82

@@ -25,6 +25,7 @@ DFLAGS += -DBASIC_SCALAR
 
 DFLAGS += -DDUST
 # DFLAGS += -DOUTFLOW_ANALYSIS
+DFLAGS += -DGLOBAL_REDUCE_DUST
 DFLAGS += -DSCALAR_FLOOR
 
 # MW_MODEL vs m82 (default)

builds/setup.summit.gcc.sh

@@ -7,5 +7,4 @@
 # module load gcc cuda fftw hdf5 python googletest/1.11.0
 module load gcc cuda/12.2.0 hdf5 python
 
-export F_OFFLOAD="-fopenmp"
 export CHOLLA_ENVSET=1

src/analysis/global_reduce_dust.cu

@@ -0,0 +1,70 @@
+#ifdef DUST
+  #ifdef GLOBAL_REDUCE_DUST
+    // STL includes
+    #include <math.h>
+    #include <stdio.h>
+
+    #include <cstdio>
+    #include <fstream>
+    #include <vector>
+
+    #include "../analysis/global_reduce_dust.h"
+    #include "../grid/grid_enum.h"
+    #include "../utils/DeviceVector.h"
+    #include "../utils/cuda_utilities.h"
+    #include "../utils/hydro_utilities.h"
+    #include "../utils/reduction_utilities.h"
+
+void Global_Reduce_Dust(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz, int n_ghost, int n_fields,
+                        Real *mass_cloud, Real *mass_dust, Real density_cloud_init)
+{
+  cuda_utilities::AutomaticLaunchParams static const launchParams(Global_Reduce_Dust_Kernel);
+
+  cuda_utilities::DeviceVector<Real> dev_mass_cloud(1);
+  cuda_utilities::DeviceVector<Real> dev_mass_dust(1);
+
+  hipLaunchKernelGGL(Global_Reduce_Dust_Kernel, launchParams.get_numBlocks(), launchParams.get_threadsPerBlock(), 0, 0,
+                     dev_conserved, nx, ny, nz, dx, dy, dz, n_ghost, n_fields, dev_mass_cloud.data(),
+                     dev_mass_dust.data(), density_cloud_init);
+  cudaDeviceSynchronize();
+
+  *mass_cloud = dev_mass_cloud[0];
+  *mass_dust  = dev_mass_dust[0];
+}
+
+__global__ void Global_Reduce_Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz,
+                                          int n_ghost, int n_fields, Real *mass_cloud, Real *mass_dust,
+                                          Real density_cloud_init)
+{
+  int xid, yid, zid, n_cells;
+  n_cells = nx * ny * nz;
+
+  Real density_gas;
+  Real mass_cloud_stride = 0.0;
+  Real density_dust;
+  Real mass_dust_stride = 0.0;
+
+  for (size_t id = threadIdx.x + blockIdx.x * blockDim.x; id < n_cells; id += blockDim.x * gridDim.x) {
+    cuda_utilities::compute3DIndices(id, nx, ny, xid, yid, zid);
+    // grid cells
+    if (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost && zid > n_ghost - 1 &&
+        zid < nz - n_ghost) {
+      density_gas  = dev_conserved[id + n_cells * grid_enum::density];
+      density_dust = dev_conserved[id + n_cells * grid_enum::dust_density];
+
+      mass_dust_stride += density_dust * dx * dy * dz;
+
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        mass_cloud_stride += density_gas * dx * dy * dz;
+      }
+    }
+  }
+
+  __syncthreads();
+
+  reduction_utilities::Grid_Reduce_Add(mass_cloud_stride, mass_cloud);
+  reduction_utilities::Grid_Reduce_Add(mass_dust_stride, mass_dust);
+}
+
+  #endif  // GLOBAL_REDUCE_DUST
+#endif    // DUST

src/analysis/global_reduce_dust.h

@@ -0,0 +1,15 @@
+#ifdef GLOBAL_REDUCE_DUST
+  #ifndef GLOBAL_REDUCE_DUST_CUDA_H
+    #define GLOBAL_REDUCE_DUST_CUDA_H
+
+    #include "../global/global.h"
+    #include "../utils/gpu.hpp"
+
+void Global_Reduce_Dust(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz, int n_ghost, int n_fields,
+                      Real *mass_cloud, Real *mass_dust, Real density_cloud_init);
+
+__global__ void Global_Reduce_Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz,
+                                        int n_ghost, int n_fields, Real *mass_cloud, Real *mass_dust,
+                                        Real density_cloud_init);
+  #endif  //  GLOBAL_REDUCE_DUST_CUDA_H
+#endif    //  GLOBAL_REDUCE_DUST

src/global/global.cpp

@@ -516,9 +516,10 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "density_wind_init") == 0) {
     parms->density_wind_init = atof(value);
 #endif
-#if defined(OUTFLOW_ANALYSIS) || defined(CLOUD_TRACKING)
+#if defined(OUTFLOW_ANALYSIS) || defined(GLOBAL_REDUCE_DUST) || defined(CLOUD_TRACKING)
   } else if (strcmp(name, "density_cloud_init") == 0) {
     parms->density_cloud_init = atof(value);
+    chprintf("Initial cloud density: %e (g/cm^3)\n", parms->density_cloud_init);
 #endif
 #ifdef SCALAR
   #ifdef DUST

src/global/global.h

@@ -312,7 +312,7 @@ struct Parameters {
 #ifdef CLOUD_TRACKING
   Real density_wind_init;
 #endif
-#if defined(CLOUD_TRACKING) || defined(OUTFLOW_ANALYSIS)
+#if defined(CLOUD_TRACKING) || defined(OUTFLOW_ANALYSIS) || defined(GLOBAL_REDUCE_DUST)
   Real density_cloud_init;
 #endif
 #ifdef TILED_INITIAL_CONDITIONS