get floor values from paramter file and remove density floor kernel

builds/make.type.dust

@@ -15,6 +15,7 @@ DFLAGS    += -DHLLC
 DFLAGS    += -DDE
 DFLAGS    += -DAVERAGE_SLOW_CELLS
 DFLAGS    += -DTEMPERATURE_FLOOR
+DFLAGS    += -DDENSITY_FLOOR
 
 DFLAGS    += -DVL
 
@@ -37,6 +38,8 @@ 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

src/global/global.cpp

@@ -428,6 +428,18 @@ 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);
+# endif
+# ifdef DENSITY_FLOOR
+  } else if (strcmp(name, "density_floor") == 0) {
+    parms->density_floor = atof(value);
+# endif
+# ifdef SCALAR_FLOOR
+  } else if (strcmp(name, "scalar_floor") == 0) {
+    parms->scalar_floor = atof(value);
+# endif
 #ifdef ANALYSIS
   } else if (strcmp(name, "analysis_scale_outputs_file") == 0) {
     strncpy(parms->analysis_scale_outputs_file, value, MAXLEN);

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
@@ -325,6 +321,15 @@ struct Parameters {
   char UVB_rates_file[MAXLEN];  // File for the UVB photoheating and
                                 // photoionization rates of HI, HeI and HeII
 #endif
+#ifdef TEMPERATURE_FLOOR
+  Real temperature_floor;
+#endif
+#ifdef DENSITY_FLOOR
+  Real density_floor;
+#endif
+#ifdef SCALAR_FLOOR
+  Real scalar_floor;
+#endif
 #ifdef ANALYSIS
   char analysis_scale_outputs_file[MAXLEN];  // File for the scale_factor output
                                              // values for cosmological

src/grid/grid3D.cpp

@@ -258,16 +258,22 @@ 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;
+#else
+  H.temperature_floor = 0.0;
+#endif
+
 #ifdef DENSITY_FLOOR
-  H.density_floor = DENS_FLOOR;
+  H.density_floor = P->density_floor;
 #else
   H.density_floor = 0.0;
 #endif
 
-#ifdef TEMPERATURE_FLOOR
-  H.temperature_floor = TEMP_FLOOR;
+#ifdef SCALAR_FLOOR
+  H.scalar_floor = P->scalar_floor;
 #else
-  H.temperature_floor = 0.0;
+  H.scalar_floor = 0.0;
 #endif
 
 #ifdef COSMOLOGY
@@ -461,12 +467,12 @@ 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,
-                         C.Grav_potential);
+                         C.Grav_potential, H.scalar_floor);
   #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);
+                             U_floor, C.Grav_potential, H.scalar_floor);
   #endif  // SIMPLE
 #endif
   } else {

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;
 

src/grid/initial_conditions.cpp

@@ -1324,7 +1324,7 @@ void Grid3D::Clouds()
   Real p_bg, p_cl;       // background and cloud pressure
   Real mu   = 0.6;       // mean atomic weight
   int N_cl  = 1;         // number of clouds
-  Real R_cl = .1;       // cloud radius in code units (kpc)
+  Real R_cl = 2.5;       // cloud radius in code units (kpc)
   Real cl_pos[N_cl][3];  // array of cloud positions
   Real r;
 
@@ -1339,26 +1339,25 @@ void Grid3D::Clouds()
 
   // single centered cloud setup
   for (int nn = 0; nn < N_cl; nn++) {
-    cl_pos[nn][0] = 0.075 * H.xdglobal;
+    cl_pos[nn][0] = 0.5 * H.xdglobal;
     cl_pos[nn][1] = 0.5 * H.ydglobal;
     cl_pos[nn][2] = 0.5 * H.zdglobal;
     printf("Cloud positions: %f %f %f\n", cl_pos[nn][0], cl_pos[nn][1], cl_pos[nn][2]);
   }
 
-  n_bg   = 1e-2;
-  n_cl   = 10;
+  n_bg   = 1.68e-4;
+  n_cl   = 5.4e-2;
   rho_bg = n_bg * mu * MP / DENSITY_UNIT;
   rho_cl = n_cl * mu * MP / DENSITY_UNIT;
-  vx_bg  = 1000*TIME_UNIT/KPC;
+  vx_bg  = 0.0;
   // vx_c  = -200*TIME_UNIT/KPC; // convert from km/s to kpc/kyr
   vx_cl = 0.0;
   vy_bg = vy_cl = 0.0;
   vz_bg = vz_cl = 0.0;
   T_bg          = 3e6;
   T_cl          = 1e4;
   p_bg          = n_bg * KB * T_bg / PRESSURE_UNIT;
-  // p_cl          = p_bg;
-  p_cl = n_cl * KB * T_cl / PRESSURE_UNIT;
+  p_cl          = p_bg;
 
   istart = H.n_ghost;
   iend   = H.nx - H.n_ghost;
@@ -1419,7 +1418,7 @@ void Grid3D::Clouds()
 #ifdef DUST
             C.host[id + H.n_cells * grid_enum::dust_density] = rho_cl * 1e-2;
 #endif  // DUST
-#endif  // SCAlAR
+#endif  // SCALAR
           }
         }
       }

src/hydro/hydro_cuda.cu

@@ -379,8 +379,7 @@ __global__ void Update_Conserved_Variables_3D(Real *dev_conserved, Real *Q_Lx, R
 
   #endif  // GRAVITY
 
-  // #if !(defined(DENSITY_FLOOR) && defined(TEMPERATURE_FLOOR))
-  #if !(defined(DENSITY_FLOOR))
+  #if !(defined(DENSITY_FLOOR) && defined(TEMPERATURE_FLOOR))
     if (dev_conserved[id] < 0.0 || dev_conserved[id] != dev_conserved[id] || dev_conserved[4 * n_cells + id] < 0.0 ||
         dev_conserved[4 * n_cells + id] != dev_conserved[4 * n_cells + id]) {
       printf("%3d %3d %3d Thread crashed in final update. %e %e %e %e %e\n", xid + x_off, yid + y_off, zid + z_off,
@@ -1140,37 +1139,6 @@ __global__ void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int
   }
 }
 
-__global__ void Apply_Density_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real density_floor)
-{
-  int id, xid, yid, zid, n_cells;
-  Real density_init;  // 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) {
-    density_init = dev_conserved[id + n_cells * grid_enum::density];
-
-    if (density_init < density_floor) {
-      printf("###Thread density change  %f -> %f \n", density_init, density_floor);
-      dev_conserved[id] = density_floor;
-      // Scale the conserved values to the new density
-      dev_conserved[id + n_cells * grid_enum::momentum_x] *= (density_floor / density_init);
-      dev_conserved[id + n_cells * grid_enum::momentum_y] *= (density_floor / density_init);
-      dev_conserved[id + n_cells * grid_enum::momentum_z] *= (density_floor / density_init);
-      dev_conserved[id + n_cells * grid_enum::Energy] *= (density_floor / density_init);
-  #ifdef DE
-      dev_conserved[id + n_cells * grid_enum::GasEnergy] *= (density_floor / density_init);
-  #endif  // DE
-    }
-  }
-}
 __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)
 {
@@ -1305,6 +1273,7 @@ __global__ void Apply_Scalar_Floor(Real *dev_conserved, int nx, int ny, int nz,
     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;
     }
   }

src/hydro/hydro_cuda.h

@@ -88,8 +88,6 @@ __global__ void Average_Slow_Cells_3D(Real *dev_conserved, int nx, int ny, int n
 __global__ void Apply_Temperature_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields,
                                         Real U_floor);
 
-__global__ void Apply_Density_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real density_floor);
-
 __global__ void Apply_Scalar_Floor(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int field_num,
                                    Real scalar_floor);
 

src/hydro/hydro_cuda_tests.cu

@@ -148,7 +148,6 @@ TEST(tMHDMhdInverseCrossingTime, CorrectInputExpectCorrectOutput)
 
 TEST(tHYDROScalarFloor, CorrectInputExpectCorrectOutput)
 {
-  // Call the function we are testing
   int num_blocks = 1;
   dim3 dim1dGrid(num_blocks, 1, 1);
   dim3 dim1dBlock(TPB, 1, 1);

src/integrators/VL_3D_cuda.cu

@@ -37,7 +37,7 @@ __global__ void Update_Conserved_Variables_3D_half(Real *dev_conserved, Real *de
 void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                           int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
                           Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
-                          Real *host_grav_potential)
+                          Real *host_grav_potential, Real scalar_floor)
 {
   // Here, *dev_conserved contains the entire
   // set of conserved variables on the grid
@@ -197,11 +197,6 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int
                      F_x, F_y, F_z, nx, ny, nz, n_ghost, dx, dy, dz, 0.5 * dt, gama, n_fields, density_floor);
   GPU_Error_Check();
 
-  #ifdef DENSITY_FLOOR
-  hipLaunchKernelGGL(Apply_Density_Floor, dim1dGrid, dim1dBlock, 0, 0, dev_conserved_half, nx, ny, nz, n_ghost,
-                     density_floor);
-  #endif  // DENSITY_FLOOR
-
   #ifdef MHD
   // Update the magnetic fields
   hipLaunchKernelGGL(mhd::Update_Magnetic_Field_3D, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, dev_conserved_half,
@@ -334,7 +329,7 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int
   #ifdef SCALAR_FLOOR
   #ifdef DUST
   hipLaunchKernelGGL(Apply_Scalar_Floor, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost,
-                     grid_enum::dust_density, 1e-10);
+                     grid_enum::dust_density, scalar_floor);
   GPU_Error_Check();
   #endif
   #endif  // SCALAR_FLOOR

src/integrators/VL_3D_cuda.h

@@ -11,7 +11,7 @@
 void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                           int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
                           Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
-                          Real *host_grav_potential);
+                          Real *host_grav_potential, Real scalar_floor);
 
 void Free_Memory_VL_3D();
 

src/integrators/simple_3D_cuda.cu

@@ -27,7 +27,7 @@
 void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                               int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
                               Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
-                              Real *host_grav_potential)
+                              Real *host_grav_potential, Real scalar_floor)
 {
   // Here, *dev_conserved contains the entire
   // set of conserved variables on the grid
@@ -189,7 +189,7 @@ void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx,
     #ifdef SCALAR_FLOOR
     #ifdef DUST
   hipLaunchKernelGGL(Apply_Scalar_Floor, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost,
-                     grid_enum::dust_density, 1e-5);
+                     grid_enum::dust_density, scalar_floor);
   CudaCheckError();
     #endif DUST
     #endif  // SCALAR_FLOOR

src/integrators/simple_3D_cuda.h

@@ -12,7 +12,7 @@
 void Simple_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                               int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
                               Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor,
-                              Real *host_grav_potential);
+                              Real *host_grav_potential, Real scalar_floor);
 
 void Free_Memory_Simple_3D();