Merge pull request #363 from helenarichie/dev-scalar-floor

Add scalar floor and refactor temperature and density floors
cholla-hydro · Feb 1, 2024 · 737d039 · 737d039
2 parents 5ff572c + cede7f3
commit 737d039
Show file tree

Hide file tree

Showing 16 changed files with 203 additions and 93 deletions.
diff --git a/builds/make.type.dust b/builds/make.type.dust
@@ -11,36 +11,36 @@ DFLAGS    += -DPRECISION=2
 DFLAGS    += -DPLMC
 DFLAGS    += -DHLLC
 
-# DFLAGS    += -DDE
+DFLAGS    += -DDE
 DFLAGS    += -DAVERAGE_SLOW_CELLS
 DFLAGS    += -DTEMPERATURE_FLOOR
+DFLAGS    += -DDENSITY_FLOOR
 
-ifeq ($(findstring cosmology,$(TYPE)),cosmology)
-DFLAGS    += -DSIMPLE
-else
 DFLAGS    += -DVL
-# DFLAGS    += -DSIMPLE
-endif
 
 # Evolve additional scalars
-DFLAGS += -DSCALAR
+DFLAGS    += -DSCALAR
+DFLAGS    += -DSCALAR_FLOOR
 
 # Define dust macro
-DFLAGS += -DDUST
+DFLAGS    += -DDUST
 
 # Apply the cooling in the GPU from precomputed tables
-DFLAGS += -DCOOLING_GPU
+DFLAGS    += -DCOOLING_GPU
+DFLAGS    += -DCLOUDY_COOLING
 
 #Measure the Timing of the different stages
-#DFLAGS += -DCPU_TIME
+#DFLAGS   += -DCPU_TIME
 
-DFLAGS += -DSLICES
-DFLAGS += -DPROJECTION
+DFLAGS    += -DSLICES
+DFLAGS    += -DPROJECTION
 
 DFLAGS    += $(OUTPUT)
 
+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)
diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu
@@ -90,8 +90,7 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
     Real const temperature = hydro_utilities::Calc_Temp_Conserved(energy, density_gas, momentum_x, momentum_y,
                                                                   momentum_z, gamma, number_density);
     #endif  // MHD
-
-  #endif  // DE
+  #endif    // DE
 
     Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) /
                   TIME_UNIT;  // sputtering timescale, kyr (sim units)
@@ -126,7 +125,6 @@ __device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real tem
   number_density /= (6e-4);  // gas number density in units of 10^-27 g/cm^3
 
   // sputtering timescale, s
-  printf("%e\n", grain_radius);
   Real tau_sp = A * (a / number_density) * (pow(temperature_0 / temperature, omega) + 1);
 
   return tau_sp;

diff --git a/src/dust/dust_cuda_tests.cpp b/src/dust/dust_cuda_tests.cpp
@@ -21,8 +21,7 @@
 
 #ifdef DUST
 
-TEST(tDUSTTestSputteringTimescale,
-     CorrectInputExpectCorrectOutput)  // test suite name, test name
+TEST(tDUSTTestSputteringTimescale, CorrectInputExpectCorrectOutput)
 {
   // Parameters
   Real YR_IN_S                     = 3.154e7;
@@ -47,8 +46,7 @@ TEST(tDUSTTestSputteringTimescale,
                        << "The ULP difference is:       " << ulps_diff << std::endl;
 }
 
-TEST(tDUSTTestSputteringGrowthRate,
-     CorrectInputExpectCorrectOutput)  // test suite name, test name
+TEST(tDUSTTestSputteringGrowthRate, CorrectInputExpectCorrectOutput)
 {
   // Parameters
   Real YR_IN_S                   = 3.154e7;

diff --git a/src/global/global.cpp b/src/global/global.cpp
@@ -436,6 +436,33 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "UVB_rates_file") == 0) {
     strncpy(parms->UVB_rates_file, value, MAXLEN);
 #endif
+#ifdef TEMPERATURE_FLOOR
+  } else if (strcmp(name, "temperature_floor") == 0) {
+    parms->temperature_floor = atof(value);
+    if (parms->temperature_floor == 0) {
+      chprintf(
+          "WARNING: temperature floor is set to its default value (zero)! It can be set to a different value in the "
+          "input parameter file.\n");
+    }
+#endif
+#ifdef DENSITY_FLOOR
+  } else if (strcmp(name, "density_floor") == 0) {
+    parms->density_floor = atof(value);
+    if (parms->density_floor == 0) {
+      chprintf(
+          "WARNING: density floor is set to its default value (zero)! It can be set to a different value in the input "
+          "parameter file.\n");
+    }
+#endif
+#ifdef SCALAR_FLOOR
+  } else if (strcmp(name, "scalar_floor") == 0) {
+    parms->scalar_floor = atof(value);
+    if (parms->scalar_floor == 0) {
+      chprintf(
+          "WARNING: scalar floor is set to its default value (zero)! It can be set to a different value in the input "
+          "parameter file.\n");
+    }
+#endif
 #ifdef ANALYSIS
   } else if (strcmp(name, "analysis_scale_outputs_file") == 0) {
     strncpy(parms->analysis_scale_outputs_file, value, MAXLEN);

diff --git a/src/global/global.h b/src/global/global.h
@@ -49,10 +49,6 @@ typedef double Real;
 
 #define LOG_FILE_NAME "run_output.log"
 
-// Conserved Floor Values
-#define TEMP_FLOOR 1e-3
-#define DENS_FLOOR 1e-5  // in code units
-
 // Parameters for Enzo dual Energy Condition
 // - Prior to GH PR #356, DE_ETA_1 nominally had a value of 0.001 in all
 //   simulations (in practice, the value of DE_ETA_1 had minimal significance
@@ -326,6 +322,9 @@ struct Parameters {
   char UVB_rates_file[MAXLEN];  // File for the UVB photoheating and
                                 // photoionization rates of HI, HeI and HeII
 #endif
+  Real temperature_floor = 0;
+  Real density_floor     = 0;
+  Real scalar_floor      = 0;
 #ifdef ANALYSIS
   char analysis_scale_outputs_file[MAXLEN];  // File for the scale_factor output
                                              // values for cosmological

diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp
@@ -258,16 +258,16 @@ void Grid3D::Initialize(struct Parameters *P)
 #endif /*ROTATED_PROJECTION*/
 
 // Values for lower limit for density and temperature
+#ifdef TEMPERATURE_FLOOR
+  H.temperature_floor = P->temperature_floor;
+#endif
+
 #ifdef DENSITY_FLOOR
-  H.density_floor = DENS_FLOOR;
-#else
-  H.density_floor = 0.0;
+  H.density_floor = P->density_floor;
 #endif
 
-#ifdef TEMPERATURE_FLOOR
-  H.temperature_floor = TEMP_FLOOR;
-#else
-  H.temperature_floor = 0.0;
+#ifdef SCALAR_FLOOR
+  H.scalar_floor = P->scalar_floor;
 #endif
 
 #ifdef COSMOLOGY
@@ -424,16 +424,6 @@ void Grid3D::Execute_Hydro_Integrator(void)
   z_off = nz_local_start;
 #endif
 
-  // Set the lower limit for density and temperature (Internal Energy)
-  Real U_floor, density_floor;
-  density_floor = H.density_floor;
-  // Minimum of internal energy from minumum of temperature
-  U_floor = H.temperature_floor * KB / (gama - 1) / MP / SP_ENERGY_UNIT;
-#ifdef COSMOLOGY
-  U_floor = H.temperature_floor / (gama - 1) / MP * KB * 1e-10;  // ( km/s )^2
-  U_floor /= Cosmo.v_0_gas * Cosmo.v_0_gas / Cosmo.current_a / Cosmo.current_a;
-#endif
-
 #ifdef CPU_TIME
   Timer.Hydro_Integrator.Start();
 #endif  // CPU_TIME
@@ -461,13 +451,13 @@ void Grid3D::Execute_Hydro_Integrator(void)
   {
 #ifdef VL
     VL_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
-                         H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, density_floor, U_floor,
+                         H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, H.density_floor,
                          C.Grav_potential);
 #endif  // VL
 #ifdef SIMPLE
     Simple_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
-                             H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, density_floor,
-                             U_floor, C.Grav_potential);
+                             H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, H.density_floor,
+                             C.Grav_potential);
 #endif  // SIMPLE
   } else {
     chprintf("Error: Grid dimensions nx: %d  ny: %d  nz: %d  not supported.\n", H.nx, H.ny, H.nz);
@@ -500,8 +490,25 @@ Real Grid3D::Update_Hydro_Grid()
 
   Execute_Hydro_Integrator();
 
-  // == Perform chemistry/cooling (there are a few different cases) ==
+#ifdef TEMPERATURE_FLOOR
+  // Set the lower limit temperature (Internal Energy)
+  Real U_floor;
+  // Minimum of internal energy from minumum of temperature
+  U_floor = H.temperature_floor * KB / (gama - 1) / MP / SP_ENERGY_UNIT;
+  #ifdef COSMOLOGY
+  U_floor = H.temperature_floor / (gama - 1) / MP * KB * 1e-10;  // ( km/s )^2
+  U_floor /= Cosmo.v_0_gas * Cosmo.v_0_gas / Cosmo.current_a / Cosmo.current_a;
+  #endif
+  Apply_Temperature_Floor(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, U_floor);
+#endif  // TEMPERATURE_FLOOR
+
+#ifdef SCALAR_FLOOR
+  #ifdef DUST
+  Apply_Scalar_Floor(C.device, H.nx, H.ny, H.nz, H.n_ghost, grid_enum::dust_density, H.scalar_floor);
+  #endif
+#endif  // SCALAR_FLOOR
 
+// == Perform chemistry/cooling (there are a few different cases) ==
 #ifdef COOLING_GPU
   #ifdef CPU_TIME
   Timer.Cooling_GPU.Start();

diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h
@@ -231,8 +231,9 @@ struct Header {
   int custom_grav;
 
   // Values for lower limit for density and temperature
-  Real density_floor;
   Real temperature_floor;
+  Real density_floor;
+  Real scalar_floor;
 
   Real Ekin_avrg;
 

diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp
@@ -1414,12 +1414,11 @@ void Grid3D::Clouds()
 #ifdef DE
             C.GasEnergy[id] = p_cl / (gama - 1.0);
 #endif  // DE
-
 #ifdef SCALAR
   #ifdef DUST
             C.host[id + H.n_cells * grid_enum::dust_density] = rho_cl * 1e-2;
   #endif  // DUST
-#endif
+#endif    // SCALAR
           }
         }
       }

diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu
@@ -1100,9 +1100,22 @@ __global__ void Sync_Energies_3D(Real *dev_conserved, int nx, int ny, int nz, in
 
 #endif  // DE
 
-#ifdef TEMPERATURE_FLOOR
-__global__ void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
-                                        Real U_floor)
+void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real U_floor)
+{
+  // set values for GPU kernels
+  int n_cells = nx * ny * nz;
+  int ngrid   = (n_cells + TPB - 1) / TPB;
+  // number of blocks per 1D grid
+  dim3 dim1dGrid(ngrid, 1, 1);
+  //  number of threads per 1D block
+  dim3 dim1dBlock(TPB, 1, 1);
+
+  hipLaunchKernelGGL(Temperature_Floor_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost,
+                     n_fields, U_floor);
+}
+
+__global__ void Temperature_Floor_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
+                                         Real U_floor)
 {
   int id, xid, yid, zid, n_cells;
   Real d, d_inv, vx, vy, vz, E, Ekin, U;
@@ -1130,15 +1143,14 @@ __global__ void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int
       dev_conserved[4 * n_cells + id] = Ekin + d * U_floor;
     }
 
-  #ifdef DE
+#ifdef DE
     U = dev_conserved[(n_fields - 1) * n_cells + id] / d;
     if (U < U_floor) {
       dev_conserved[(n_fields - 1) * n_cells + id] = d * U_floor;
     }
-  #endif
+#endif
   }
 }
-#endif  // TEMPERATURE_FLOOR
 
 __device__ Real Average_Cell_Single_Field(int field_indx, int i, int j, int k, int nx, int ny, int nz, int ncells,
                                           Real *conserved)
@@ -1254,3 +1266,42 @@ __device__ void Average_Cell_All_Fields(int i, int j, int k, int nx, int ny, int
 
   printf("%3d %3d %3d FC: d: %e  E:%e  P:%e  vx:%e  vy:%e  vz:%e\n", i, j, k, d, E, P, vx_av, vy_av, vz_av);
 }
+
+void Apply_Scalar_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int field_num, Real scalar_floor)
+{
+  // set values for GPU kernels
+  int n_cells = nx * ny * nz;
+  int ngrid   = (n_cells + TPB - 1) / TPB;
+  // number of blocks per 1D grid
+  dim3 dim1dGrid(ngrid, 1, 1);
+  //  number of threads per 1D block
+  dim3 dim1dBlock(TPB, 1, 1);
+
+  hipLaunchKernelGGL(Scalar_Floor_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, field_num,
+                     scalar_floor);
+}
+
+__global__ void Scalar_Floor_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int field_num,
+                                    Real scalar_floor)
+{
+  int id, xid, yid, zid, n_cells;
+  Real scalar;  // variable to store the value of the scalar before a floor is applied
+  n_cells = nx * ny * nz;
+
+  // get a global thread ID
+  id  = threadIdx.x + blockIdx.x * blockDim.x;
+  zid = id / (nx * ny);
+  yid = (id - zid * nx * ny) / nx;
+  xid = id - zid * nx * ny - yid * nx;
+
+  // threads corresponding to real cells do the calculation
+  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) {
+    scalar = dev_conserved[id + n_cells * field_num];
+
+    if (scalar < scalar_floor) {
+      // printf("###Thread scalar change  %f -> %f \n", scalar, scalar_floor);
+      dev_conserved[id + n_cells * field_num] = scalar_floor;
+    }
+  }
+}
diff --git a/src/hydro/hydro_cuda.h b/src/hydro/hydro_cuda.h
@@ -84,10 +84,15 @@ __global__ void Average_Slow_Cells_3D(Real *dev_conserved, int nx, int ny, int n
                                       Real dy, Real dz, Real gamma, Real max_dti_slow);
 #endif
 
-#ifdef TEMPERATURE_FLOOR
-__global__ void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
-                                        Real U_floor);
-#endif
+void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real U_floor);
+
+__global__ void Temperature_Floor_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
+                                         Real U_floor);
+
+void Apply_Scalar_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int field_num, Real scalar_floor);
+
+__global__ void Scalar_Floor_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int field_num,
+                                    Real scalar_floor);
 
 __global__ void Partial_Update_Advected_Internal_Energy_1D(Real *dev_conserved, Real *Q_Lx, Real *Q_Rx, int nx,
                                                            int n_ghost, Real dx, Real dt, Real gamma, int n_fields);