resolve conflicts for destroyed cloud simulation

helenarichie · Jun 13, 2024 · 48817d7 · 48817d7
2 parents 24ffdab + 631c354
commit 48817d7
Show file tree

Hide file tree

Showing 31 changed files with 256 additions and 203 deletions.
diff --git a/builds/make.type.dust b/builds/make.type.dust
@@ -11,16 +11,16 @@ DFLAGS    += -DPRECISION=2
 DFLAGS    += -DPPMP
 DFLAGS    += -DHLLC
 
-DFLAGS    += -DDE
-DFLAGS    += -DAVERAGE_SLOW_CELLS
+# DFLAGS    += -DDE
+# DFLAGS    += -DAVERAGE_SLOW_CELLS
 DFLAGS    += -DTEMPERATURE_FLOOR
-DFLAGS    += -DDENSITY_FLOOR
+# DFLAGS    += -DDENSITY_FLOOR
 
 DFLAGS    += -DVL
 
 # Evolve additional scalars
-DFLAGS    += -DSCALAR
-DFLAGS    += -DSCALAR_FLOOR
+# DFLAGS    += -DSCALAR
+# DFLAGS    += -DSCALAR_FLOOR
 
 # Define dust macro
 DFLAGS    += -DDUST
@@ -35,8 +35,8 @@ DFLAGS    += -DCOOLING_GPU
 #Measure the Timing of the different stages
 #DFLAGS   += -DCPU_TIME
 
-DFLAGS    += -DSLICES
-DFLAGS    += -DPROJECTION
+# DFLAGS    += -DSLICES
+# DFLAGS    += -DPROJECTION
 
 DFLAGS    += $(OUTPUT)
 
@@ -46,4 +46,4 @@ DFLAGS    += $(OUTPUT)
 #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)
diff --git a/builds/setup.summit.gcc.sh b/builds/setup.summit.gcc.sh
@@ -4,7 +4,8 @@
 #   source ./setup.summit.gcc.sh
 
 #module load gcc/10.2.0 cuda/11.4.0 fftw hdf5 python
-module load gcc cuda fftw hdf5 python googletest/1.11.0
+#module load gcc cuda fftw hdf5 python googletest/1.11.0
+# module load gcc cuda/12.2.0 fftw hdf5 python
+module load gcc cuda/12.2.0 hdf5 python ums/default ums-ompix/default openmpi/5.0.3
 
-export F_OFFLOAD="-fopenmp"
 export CHOLLA_ENVSET=1
diff --git a/cloud-wind-surv.txt b/cloud-wind-surv.txt
@@ -0,0 +1,44 @@
+################################################
+# Parameters required for all problems
+################################################
+# number of grid cells in the x dimension
+nx=1024
+# number of grid cells in the y dimension
+ny=128
+# number of grid cells in the z dimension
+nz=128
+# x direction lower domain boundary
+xmin=0.0
+# x direction domain length
+xlen=12.8
+# y direction lower domain boundary
+ymin=0.0
+# y direction domain length
+ylen=1.6
+# z direction lower domain boundary
+zmin=0
+# y direction domain length
+zlen=1.6
+# final output time
+tout=120000
+# time interval for output
+outstep=50
+n_hydro=5
+# ratio of specific heats
+gamma=1.66666667
+# name of initial conditions
+init=Clouds
+# type of boundary condition, options include 1 (periodic), 2 (reflective), 3 (transmissive), 4 (custom)
+xl_bcnd=4
+xu_bcnd=3
+yl_bcnd=3
+yu_bcnd=3
+zl_bcnd=3
+zu_bcnd=3
+# path to output directory
+outdir=./hdf5/raw/
+custom_bcnd=wind
+temperature_floor=10
+scalar_floor=1e-10
+density_cloud_init=1e-23
+grain_radius=1
diff --git a/src/analysis/feedback_analysis.cpp b/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++) {

diff --git a/src/analysis/outflow_analysis.cu b/src/analysis/outflow_analysis.cu
@@ -1,54 +1,55 @@
 #ifdef DUST
-#ifdef OUTFLOW_ANALYSIS
-  // STL includes
-  #include <stdio.h>
-  #include <math.h> 
-
-  #include <cstdio>
-  #include <fstream>
-  #include <vector>
-
-  #include "../grid/grid_enum.h"
-  #include "../utils/DeviceVector.h"
-  #include "../utils/cuda_utilities.h"
-  #include "../utils/reduction_utilities.h"
-  #include "../analysis/outflow_analysis.h"
-  #include "../utils/hydro_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)
+  #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.numBlocks, launchParams.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);
+  hipLaunchKernelGGL(Outflow_Analysis_Kernel, launchParams.numBlocks, launchParams.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];
+  *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)
+__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;
+  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_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;
@@ -58,7 +59,7 @@ __global__ void Outflow_Analysis_Kernel(Real *dev_conserved, int nx, int ny, int
       }
     }
   }
-  
+
   __syncthreads();
 
   reduction_utilities::Grid_Reduce_Add(mass_cloud_stride, mass_cloud);
@@ -236,5 +237,5 @@ __global__ void Outflow_Analysis_Kernel(Real *dev_conserved, int nx, int ny, int
   */
 }
 
-#endif  // OUTFLOW_ANALYSIS
-#endif  // DUST
+  #endif  // OUTFLOW_ANALYSIS
+#endif    // DUST
diff --git a/src/analysis/outflow_analysis.h b/src/analysis/outflow_analysis.h
@@ -5,12 +5,11 @@
     #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);
 
-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);
+__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
-
diff --git a/src/chemistry_gpu/chemistry_functions.cpp b/src/chemistry_gpu/chemistry_functions.cpp
@@ -228,7 +228,7 @@ void Grid3D::Update_Chemistry()
   #ifdef COSMOLOGY
   Chem.H.current_z = Cosmo.current_z;
   #else
-  Chem.H.current_z = 0;
+  Chem.H.current_z          = 0;
   #endif
 
   Do_Chemistry_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, Chem.H);

diff --git a/src/cloud_tracking/cloud_tracking.cu b/src/cloud_tracking/cloud_tracking.cu
@@ -78,7 +78,7 @@ __global__ void Cloud_Reduction_Kernel(Real *dev_conserved, int nx, int ny, int
       velocity_x = dev_conserved[id + n_cells * grid_enum::momentum_x] / density;
       mass       = density * dx * dy * dz;
       // if ((density * DENSITY_UNIT) >= (pow(density_cloud_init * density_wind_init, 0.5))) {
-      if ((density * DENSITY_UNIT) >= (density_cloud_init / pow(1e3, 1/2))) {
+      if ((density * DENSITY_UNIT) >= (density_cloud_init / pow(1e3, 1 / 2))) {
         mass_stride += mass;
         // (Shin et al. (2008) eq. 9)
         integrand_stride += velocity_x * density * dx * dy * dz;

diff --git a/src/cooling/cooling_cuda.cu b/src/cooling/cooling_cuda.cu
@@ -107,60 +107,61 @@ __global__ void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int
   #ifdef DE
     T_init = d * ge * (gamma - 1.0) * PRESSURE_UNIT / (n * KB);
   #endif
-    // calculate cooling rate per volume
-
-    T = T_init;
+    if (T_init < 1e6) {
+      // calculate cooling rate per volume
+      T = T_init;
   // call the cooling function
   #ifdef CLOUDY_COOL
-    cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
+      cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
   #else
-    cool = CIE_cool(n, T);
+      cool = CIE_cool(n, T);
   #endif
 
-    // calculate change in temperature given dt
-    del_T = cool * dt * TIME_UNIT * (gamma - 1.0) / (n * KB);
-
-    // limit change in temperature to 1%
-    while (del_T / T > 0.01) {
-      // what dt gives del_T = 0.01*T?
-      dt_sub = 0.01 * T * n * KB / (cool * TIME_UNIT * (gamma - 1.0));
-      // apply that dt
-      T -= cool * dt_sub * TIME_UNIT * (gamma - 1.0) / (n * KB);
-      // how much time is left from the original timestep?
-      dt -= dt_sub;
+      // calculate change in temperature given dt
+      del_T = cool * dt * TIME_UNIT * (gamma - 1.0) / (n * KB);
+
+      // limit change in temperature to 1%
+      while (del_T / T > 0.01) {
+        // what dt gives del_T = 0.01*T?
+        dt_sub = 0.01 * T * n * KB / (cool * TIME_UNIT * (gamma - 1.0));
+        // apply that dt
+        T -= cool * dt_sub * TIME_UNIT * (gamma - 1.0) / (n * KB);
+        // how much time is left from the original timestep?
+        dt -= dt_sub;
   // calculate cooling again
   #ifdef CLOUDY_COOL
-      cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
+        cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
   #else
-      cool = CIE_cool(n, T);
+        cool = CIE_cool(n, T);
   #endif
-      // calculate new change in temperature
-      del_T = cool * dt * TIME_UNIT * (gamma - 1.0) / (n * KB);
-    }
+        // calculate new change in temperature
+        del_T = cool * dt * TIME_UNIT * (gamma - 1.0) / (n * KB);
+      }
 
-    // calculate final temperature
-    T -= del_T;
+      // calculate final temperature
+      T -= del_T;
 
-    // adjust value of energy based on total change in temperature
-    del_T = T_init - T;  // total change in T
-    E -= n * KB * del_T / ((gamma - 1.0) * ENERGY_UNIT);
+      // adjust value of energy based on total change in temperature
+      del_T = T_init - T;  // total change in T
+      E -= n * KB * del_T / ((gamma - 1.0) * ENERGY_UNIT);
   #ifdef DE
-    ge -= KB * del_T / (mu * MP * (gamma - 1.0) * SP_ENERGY_UNIT);
+      ge -= KB * del_T / (mu * MP * (gamma - 1.0) * SP_ENERGY_UNIT);
   #endif
 
   // calculate cooling rate for new T
   #ifdef CLOUDY_COOL
-    cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
+      cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
   #else
-    cool = CIE_cool(n, T);
-    // printf("%d %d %d %e %e %e\n", xid, yid, zid, n, T, cool);
+      cool = CIE_cool(n, T);
+  // printf("%d %d %d %e %e %e\n", xid, yid, zid, n, T, cool);
   #endif
 
-    // and send back from kernel
-    dev_conserved[4 * n_cells + id] = E;
+      // and send back from kernel
+      dev_conserved[4 * n_cells + id] = E;
   #ifdef DE
-    dev_conserved[(n_fields - 1) * n_cells + id] = d * ge;
+      dev_conserved[(n_fields - 1) * n_cells + id] = d * ge;
   #endif
+    }
   }
 }
 

diff --git a/src/cooling_grackle/cool_grackle.cpp b/src/cooling_grackle/cool_grackle.cpp
@@ -89,7 +89,7 @@ void Cool_GK::Initialize(struct Parameters *P, Cosmology &Cosmo)
   data->metal_cooling = 1;  // metal cooling off
   #else
   chprintf("WARNING: Metal Cooling is Off. \n");
-  data->metal_cooling = 0;  // metal cooling off
+  data->metal_cooling       = 0;  // metal cooling off
   #endif
 
   #ifdef PARALLEL_OMP