cholla-hydro · alwinm · Feb 6, 2023 · Feb 9, 2023 · Apr 6, 2023 · Apr 11, 2023
diff --git a/Makefile b/Makefile
@@ -9,16 +9,10 @@ include builds/make.type.$(TYPE)
 # CUDA_ARCH defaults to sm_70 if not set in make.host
 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/model src/mpi src/old_cholla src/particles src/reconstruction \
-            src/riemann_solvers src/system_tests src/utils src/dust
-
 SUFFIX ?= .$(TYPE).$(MACHINE)
 
-CPPFILES := $(foreach DIR,$(DIRS),$(wildcard $(DIR)/*.cpp))
-GPUFILES := $(foreach DIR,$(DIRS),$(wildcard $(DIR)/*.cu))
+CPPFILES := $(shell find src/ -type f -name '*.cpp')
+GPUFILES := $(shell find src/ -type f -name '*.cu')
 
 # Build a list of all potential object files so cleaning works properly
 CLEAN_OBJS := $(subst .cpp,.o,$(CPPFILES)) \
@@ -97,7 +91,7 @@ ifeq ($(findstring -DPARIS,$(DFLAGS)),-DPARIS)
   endif
 endif
 
-ifeq ($(findstring -DSUPERNOVA,$(DFLAGS)),-DSUPERNOVA)
+ifeq ($(findstring -DFEEDBACK,$(DFLAGS)),-DFEEDBACK)
     ifdef HIPCONFIG
 	CXXFLAGS += -I$(ROCM_PATH)/include/hiprand -I$(ROCM_PATH)/hiprand/include
 	GPUFLAGS += -I$(ROCM_PATH)/include/hiprand -I$(ROCM_PATH)/hiprand/include

diff --git a/builds/make.type.disk b/builds/make.type.disk
@@ -1,12 +1,12 @@
 MPI_GPU = -DMPI_GPU
 DFLAGS += -DPARTICLES
-#DFLAGS += -DPARTICLES_CPU
 DFLAGS += -DPARTICLES_GPU
-#DFLAGS += -DONLY_PARTICLES
 DFLAGS += -DPARTICLE_IDS
-#DFLAGS += -DSINGLE_PARTICLE_MASS
 DFLAGS += -DPARTICLE_AGE
-DFLAGS += -DSUPERNOVA  #this flag requires PARTICLE_AGE, PARTICLE_IDS
+DFLAGS += -DFEEDBACK  #this flag requires PARTICLE_AGE, PARTICLE_IDS
+#DFLAGS += -DONLY_RESOLVED
+#DFLAGS += -DNO_WIND_FEEDBACK
+#DFLAGS += -DNO_SN_FEEDBACK
 DFLAGS += -DANALYSIS
 #DFLAGS += -DPARTICLES_KDK
 
@@ -29,16 +29,13 @@ DFLAGS    += -DPPMC
 DFLAGS    += -DHLLC
 DFLAGS    += -DVL
 
-DFLAGS += -DDISK_ICS
-
 DFLAGS    += -DDENSITY_FLOOR
 DFLAGS    += -DTEMPERATURE_FLOOR
 DFLAGS    += -DCOOLING_GPU
 #DFLAGS    += -DCLOUDY_COOL
 DFLAGS    += -DDE
-DFLAGS    += -DCPU_TIME
+#DFLAGS    += -DCPU_TIME
 DFLAGS    += -DAVERAGE_SLOW_CELLS
-DFLAGS    += -DHYDRO_GPU
 
 OUTPUT    ?=  -DOUTPUT -DHDF5 -DSLICES -DPROJECTION
 DFLAGS    += $(OUTPUT)
@@ -48,4 +45,4 @@ DFLAGS    += $(MPI_GPU)
 DFLAGS += -DPARALLEL_OMP
 DFLAGS += -DN_OMP_THREADS=$(OMP_NUM_THREADS)
 
-#DFLAGS += -DCUDA_ERROR_CHECK
+DFLAGS += -DCUDA_ERROR_CHECK
diff --git a/src/analysis/feedback_analysis.cpp b/src/analysis/feedback_analysis.cpp
@@ -9,8 +9,13 @@
 
 #define VRMS_CUTOFF_DENSITY (0.01 * 0.6 * MP / DENSITY_UNIT)
 
-FeedbackAnalysis::FeedbackAnalysis(Grid3D& G)
+FeedbackAnalysis::FeedbackAnalysis(Grid3D& G, struct parameters* P)
 {
+  // set distance limits beyond which contributions to the V_rms don't
+  // make sense, such as too close to the simulation volume edge or
+  // too high above the disk.
+  r_max = P->xlen / 2 * 0.975;
+  z_max = 0.15;  // 150 pc above/below the disk plane
   // allocate arrays
   h_circ_vel_x = (Real*)malloc(G.H.n_cells * sizeof(Real));
   h_circ_vel_y = (Real*)malloc(G.H.n_cells * sizeof(Real));

diff --git a/src/analysis/feedback_analysis.h b/src/analysis/feedback_analysis.h
@@ -8,6 +8,7 @@
 class FeedbackAnalysis
 {
   Real *h_circ_vel_x, *h_circ_vel_y;
+  Real r_max, z_max;
 
 #ifdef PARTICLES_GPU
   Real *d_circ_vel_x, *d_circ_vel_y;
@@ -21,8 +22,10 @@ class FeedbackAnalysis
   Real totalEnergy{0};
   Real totalMomentum{0};
   Real totalUnresEnergy{0};
+  Real totalWindMomentum{0};
+  Real totalWindEnergy{0};
 
-  FeedbackAnalysis(Grid3D& G);
+  FeedbackAnalysis(Grid3D& G, struct parameters* P);
   ~FeedbackAnalysis();
 
   void Compute_Gas_Velocity_Dispersion(Grid3D& G);

diff --git a/src/analysis/feedback_analysis_gpu.cu b/src/analysis/feedback_analysis_gpu.cu
@@ -6,7 +6,6 @@
 #include "feedback_analysis.h"
 #ifdef PARTICLES_GPU
 
-  #define MU 0.6
   // in cgs, this is 0.01 cm^{-3}
   #define MIN_DENSITY  (0.01 * MP * MU * LENGTH_UNIT * LENGTH_UNIT * LENGTH_UNIT / MASS_UNIT)  // 148279.7
   #define TPB_ANALYSIS 1024
@@ -33,7 +32,9 @@ __device__ void warpReduce(volatile Real *buff, size_t tid)
   }
 }
 
-void __global__ Reduce_Tubulence_kernel(int nx, int ny, int nz, int n_ghost, Real *density, Real *momentum_x,
+void __global__ Reduce_Tubulence_kernel(int nx, int ny, int nz, int n_ghost, Real xbound, Real ybound, Real zbound,
+                                        Real dx, Real dy, Real dz, Real nx_local_start, Real ny_local_start,
+                                        Real nz_local_start, Real r_max, Real z_max, Real *density, Real *momentum_x,
                                         Real *momentum_y, Real *momentum_z, Real *circ_vel_x, Real *circ_vel_y,
                                         Real *partial_mass, Real *partial_vel)
 {
@@ -50,15 +51,24 @@ void __global__ Reduce_Tubulence_kernel(int nx, int ny, int nz, int n_ghost, Rea
   s_mass[tid] = 0;
   s_vel[tid]  = 0;
   Real vx, vy, vz;
+  Real x, y, z, r;
+
   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[id] > MIN_DENSITY) {
-    s_mass[tid] = density[id];
-    vx          = momentum_x[id] / density[id];
-    vy          = momentum_y[id] / density[id];
-    vz          = momentum_z[id] / density[id];
-    s_vel[tid] =
-        ((vx - circ_vel_x[id]) * (vx - circ_vel_x[id]) + (vy - circ_vel_y[id]) * (vy - circ_vel_y[id]) + (vz * vz)) *
-        density[id];
+    x = xbound + (nx_local_start + xid - n_ghost + 0.5) * dx;
+    y = ybound + (ny_local_start + yid - n_ghost + 0.5) * dy;
+    z = zbound + (nz_local_start + zid - n_ghost + 0.5) * dz;
+    r = sqrt(x * x + y * y + z * z);
+
+    if (r <= r_max && abs(z) <= z_max) {
+      s_mass[tid] = density[id];
+      vx          = momentum_x[id] / density[id];
+      vy          = momentum_y[id] / density[id];
+      vz          = momentum_z[id] / density[id];
+      s_vel[tid] =
+          ((vx - circ_vel_x[id]) * (vx - circ_vel_x[id]) + (vy - circ_vel_y[id]) * (vy - circ_vel_y[id]) + (vz * vz)) *
+          density[id];
+    }
   }
   __syncthreads();
 
@@ -153,9 +163,17 @@ void FeedbackAnalysis::Compute_Gas_Velocity_Dispersion_GPU(Grid3D &G)
   Real total_mass = 0;
   Real total_vel  = 0;
 
+  Real nx_loc_strt = 0, ny_loc_strt = 0, nz_loc_strt = 0;
+  #ifdef MPI_CHOLLA
+  nx_loc_strt = nx_local_start;
+  ny_loc_strt = ny_local_start;
+  nz_loc_strt = nz_local_start;
+  #endif
+
   hipLaunchKernelGGL(Reduce_Tubulence_kernel, ngrid, TPB_ANALYSIS, 0, 0, G.H.nx, G.H.ny, G.H.nz, G.H.n_ghost,
-                     G.C.d_density, G.C.d_momentum_x, G.C.d_momentum_y, G.C.d_momentum_z, d_circ_vel_x, d_circ_vel_y,
-                     d_partial_mass, d_partial_vel);
+                     G.H.xbound, G.H.ybound, G.H.zbound, G.H.dx, G.H.dy, G.H.dz, nx_loc_strt, ny_loc_strt, nz_loc_strt,
+                     r_max, z_max, G.C.d_density, G.C.d_momentum_x, G.C.d_momentum_y, G.C.d_momentum_z, d_circ_vel_x,
+                     d_circ_vel_y, d_partial_mass, d_partial_vel);
 
   size_t n         = ngrid;
   Real *mass_input = d_partial_mass;