save m82 dust setup work in progress

.gitmodules

@@ -1,3 +0,0 @@
-[submodule "cholla-tests-data"]
-	path = cholla-tests-data
-	url = https://github.com/cholla-hydro/cholla-tests-data.git

Makefile

@@ -11,9 +11,9 @@ CUDA_ARCH ?= sm_70
 
 DIRS     := src src/analysis src/chemistry_gpu src/cooling src/cooling_grackle src/cosmology \
             src/cpu src/global src/gravity src/gravity/paris src/grid src/hydro \
-            src/integrators src/io src/main.cpp src/main_tests.cpp src/mhd\
+            src/integrators src/io src/main.cpp src/main_tests.cpp src/mhd \
             src/model src/mpi src/old_cholla src/particles src/reconstruction \
-            src/riemann_solvers src/system_tests src/utils src/dust src/cluster \
+            src/riemann_solvers src/system_tests src/utils src/dust src/cluster src/cloud_tracking \
 
 SUFFIX ?= .$(TYPE).$(MACHINE)
 

builds/make.type.dust

@@ -8,7 +8,7 @@ MPI_GPU   ?=
 
 DFLAGS    += -DMPI_CHOLLA
 DFLAGS    += -DPRECISION=2
-DFLAGS    += -DPLMC
+DFLAGS    += -DPPMP
 DFLAGS    += -DHLLC
 
 DFLAGS    += -DDE
@@ -24,10 +24,11 @@ DFLAGS    += -DSCALAR_FLOOR
 
 # Define dust macro
 DFLAGS    += -DDUST
+DFLAGS    += -DOUTFLOW_ANALYSIS
+# DFLAGS    += -DCLOUD_TRACKING
 
 # Apply the cooling in the GPU from precomputed tables
 DFLAGS    += -DCOOLING_GPU
-DFLAGS    += -DCLOUDY_COOLING
 
 #Measure the Timing of the different stages
 #DFLAGS   += -DCPU_TIME
@@ -37,10 +38,10 @@ DFLAGS    += -DPROJECTION
 
 DFLAGS    += $(OUTPUT)
 
-DFLAGS    += -DOUTPUT_ALWAYS
+# DFLAGS    += -DOUTPUT_ALWAYS
 
 #Select if the Hydro Conserved data will reside in the GPU
 #and the MPI transfers are done from the GPU
 #If not specified, MPI_GPU is off by default
 #This is set in the system make.host file
-DFLAGS    += $(MPI_GPU)
+DFLAGS    += $(MPI_GPU)

builds/make.type.hydro

@@ -13,19 +13,28 @@ DFLAGS    += -DVL
 endif
 
 # Apply a density and temperature floor
-DFLAGS    += -DDENSITY_FLOOR
-DFLAGS    += -DTEMPERATURE_FLOOR
+#DFLAGS    += -DDENSITY_FLOOR
+#DFLAGS    += -DTEMPERATURE_FLOOR
 
 # Solve the Gas Internal Energy usisng a Dual Energy Formalism
-#DFLAGS    += -DDE
+DFLAGS    += -DDE
+
+DFLAGS += -DSLICES
 
 # Apply cooling on the GPU from precomputed tables
 #DFLAGS    += -DCOOLING_GPU
+DFLAGS    += -DAVERAGE_SLOW_CELLS
+
+# Apply reference frame shift to track cloud
+DFLAGS += -DCLOUD_TRACKING
 
 # Measure the Timing of the different stages
 #DFLAGS    += -DCPU_TIME
 
 # Select output format
-# Can also add -DSLICES and -DPROJECTIONS
-OUTPUT    ?=  -DOUTPUT -DHDF5
+# Can also add -DSLICES and -DPROJECTION
+OUTPUT    ?=  -DOUTPUT -DHDF5 
 DFLAGS    += $(OUTPUT)
+
+DFLAGS    += -DSLICES
+DFLAGS    += -DPROJECTION

builds/make.type.m82

@@ -2,7 +2,9 @@
 
 #-- separated output flag so that it can be overriden in target-specific 
 #   for make check
-OUTPUT    ?=  -DOUTPUT -DHDF5 -DSLICES -DPROJECTION -DROTATED_PROJECTION
+# OUTPUT    ?=  -DOUTPUT -DHDF5 -DSLICES -DPROJECTION -DROTATED_PROJECTION
+
+OUTPUT    ?=  -DOUTPUT -DHDF5 -DSLICES
 
 MPI_GPU   ?=
 
@@ -19,10 +21,14 @@ DFLAGS    += -DSTATIC_GRAV
 
 # EVOLVE SCALAR (INJECTED BY SUPERNOVA)
 DFLAGS += -DSCALAR
-DFLAGS += -DBASIC_SCALAR
+# DFLAGS += -DBASIC_SCALAR
+
+DFLAGS += -DDUST
+# DFLAGS += -DOUTFLOW_ANALYSIS
+DFLAGS += -DSCALAR_FLOOR
 
 # MW_MODEL vs m82 (default)
-DFLAGS += -DMW_MODEL
+# DFLAGS += -DMW_MODEL
 
 ifeq ($(findstring cosmology,$(TYPE)),cosmology)
 DFLAGS    += -DSIMPLE

builds/make.type.particles

@@ -11,11 +11,12 @@ DFLAGS += -DPARTICLES_GPU
 
 
 #Solve only an N-Body Simulation (No Hydro, DM Only simulation)
-# DFLAGS += -DONLY_PARTICLES
+DFLAGS += -DONLY_PARTICLES
 
 
 # Track Particles IDs and write them to the output files
 DFLAGS += -DPARTICLE_IDS
+DFLAGS += -DCLOUD_TRACKING
 
 
 # Track Particle Ages ( Stellar Populations )

src/analysis/feedback_analysis.cpp

@@ -87,7 +87,7 @@ void FeedbackAnalysis::Compute_Gas_Velocity_Dispersion(Grid3D& G)
   #ifdef MPI_CHOLLA
   MPI_Allreduce(&partial_mass, &total_mass, 1, MPI_CHREAL, MPI_SUM, world);
   #else
-  total_mass = partial_mass;
+  total_mass         = partial_mass;
   #endif
 
   for (k = G.H.n_ghost; k < G.H.nz - G.H.n_ghost; k++) {

src/analysis/outflow_analysis.cu

@@ -0,0 +1,241 @@
+#ifdef DUST
+  #ifdef OUTFLOW_ANALYSIS
+    // STL includes
+    #include <math.h>
+    #include <stdio.h>
+
+    #include <cstdio>
+    #include <fstream>
+    #include <vector>
+
+    #include "../analysis/outflow_analysis.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 Outflow_Analysis(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(Outflow_Analysis_Kernel);
+
+  cuda_utilities::DeviceVector<Real> dev_mass_cloud(1);
+  cuda_utilities::DeviceVector<Real> dev_mass_dust(1);
+
+  hipLaunchKernelGGL(Outflow_Analysis_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 Outflow_Analysis_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 / 30)) {
+        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);
+
+  /*
+  // int id, xid, yid, zid, gid;
+
+  // calculate ghost cell ID and i,j,k in GPU grid
+  int id = threadIdx.x + blockIdx.x * blockDim.x;
+
+  int isize, jsize;
+
+  // +x boundary first
+  isize = n_ghost;
+  jsize = ny;
+
+  // not true i,j,k but relative i,j,k in the GPU grid
+  zid = id / (isize * jsize);
+  yid = (id - zid * isize * jsize) / isize;
+  xid = id - zid * isize * jsize - yid * isize;
+
+  // map thread id to ghost cell id
+  xid += nx - n_ghost;  // +x boundary
+  int gid = xid + yid * nx + zid * nx * ny;
+
+  // boundary cells
+  // -x boundary
+  if (xid == n_ghost && xid < nx && yid < ny && zid < nz) {
+    velocity_x = dev_conserved[gid + n_cells * grid_enum::momentum_x] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_x < 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        printf("-x %e/n", density_gas * DENSITY_UNIT);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_x) * dx * dx;
+      }
+      #ifdef DUST
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_x) * dx * dx;
+      #endif  // DUST
+    }
+  }
+  // +x boundary
+  if (xid == nx - n_ghost && xid < nx && yid < ny && zid < nz) {
+    velocity_x = dev_conserved[gid + n_cells * grid_enum::momentum_x] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_x > 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        printf("+x %e/n", density_gas * DENSITY_UNIT);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_x) * dx * dx;
+      }
+      #ifdef DUST
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_x) * dx * dx;
+      #endif  // DUST
+    }
+  }
+
+  // +y boundary next
+  isize = nx;
+  jsize = n_ghost;
+  // ksize = nz;
+
+  // not true i,j,k but relative i,j,k
+  zid = id / (isize * jsize);
+  yid = (id - zid * isize * jsize) / isize;
+  xid = id - zid * isize * jsize - yid * isize;
+
+  // map thread id to ghost cell id
+  yid += ny - n_ghost;
+  gid = xid + yid * nx + zid * nx * ny;
+
+  // -y boundary
+  if (xid < nx && yid == n_ghost && yid < ny && zid < nz) {
+    velocity_y = dev_conserved[gid + n_cells * grid_enum::momentum_y] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_y < 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        // printf("-y cloud %e/n", density_gas * DENSITY_UNIT);
+        printf("yid: %d\n", yid);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_y) * dx * dx;
+      }
+      #ifdef DUST
+      // printf("-y dust %e/n", density_dust * DENSITY_UNIT);
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_y) * dx * dx;
+      #endif  // DUST
+    }
+  }
+  // +y boundary
+  if (xid < nx && yid == ny - n_ghost && yid < ny && zid < nz) {
+    velocity_y = dev_conserved[gid + n_cells * grid_enum::momentum_y] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_y > 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        // printf("+y %e/n", density_gas * DENSITY_UNIT);
+        printf("yid: %d\n", yid);
+        printf("yid: %d\n", ny_real);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_y) * dx * dx;
+      }
+      #ifdef DUST
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_y) * dx * dx;
+      #endif  // DUST
+    }
+  }
+
+  isize = nx;
+  jsize = ny;
+  // ksize = n_ghost;
+
+  // not true i,j,k but relative i,j,k
+  zid = id / (isize * jsize);
+  yid = (id - zid * isize * jsize) / isize;
+  xid = id - zid * isize * jsize - yid * isize;
+
+  // map thread id to ghost cell id
+  zid += nz - n_ghost;
+  gid = xid + yid * nx + zid * nx * ny;
+
+  // -z boundary
+  if (xid < nx && yid < ny && zid == n_ghost && zid < nz) {
+    velocity_z = dev_conserved[gid + n_cells * grid_enum::momentum_z] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_z < 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        printf("-z %e/n", density_gas * DENSITY_UNIT);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_z) * dx * dx;
+      }
+      #ifdef DUST
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_z) * dx * dx;
+      #endif  // DUST
+    }
+  }
+  // +z bondary
+  if (xid < nx && yid < ny && zid == nz - n_ghost && zid < nz) {
+    velocity_z = dev_conserved[gid + n_cells * grid_enum::momentum_z] / density_gas;
+    density_gas = dev_conserved[gid + n_cells * grid_enum::density];
+    #ifdef DUST
+    density_dust = dev_conserved[gid + n_cells * grid_enum::dust_density];
+    #endif
+    if (velocity_z > 0) {
+      if ((density_gas * DENSITY_UNIT) >= (density_cloud_init / 3)) {
+        printf("+z %e/n", density_gas * DENSITY_UNIT);
+        mass_cloud_bndry_stride += density_gas * dx * dy * dz;
+        rate_cloud_stride += density_gas * std::abs(velocity_z) * dx * dx;
+      }
+      #ifdef DUST
+      mass_dust_bndry_stride += density_dust * dx * dy * dz;
+      rate_dust_stride += density_dust * std::abs(velocity_z) * dx * dx;
+      #endif  // DUST
+    }
+  }
+  */
+}
+
+  #endif  // OUTFLOW_ANALYSIS
+#endif    // DUST

src/analysis/outflow_analysis.h

@@ -0,0 +1,15 @@
+#ifdef OUTFLOW_ANALYSIS
+  #ifndef OUTFLOW_ANALYSIS_CUDA_H
+    #define OUTFLOW_ANALYSIS_CUDA_H
+
+    #include "../global/global.h"
+    #include "../utils/gpu.hpp"
+
+void Outflow_Analysis(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 Outflow_Analysis_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  //  OUTFLOW_ANALYSIS_H
+#endif    //  OUTFLOW_ANALYSIS