From 094f1e9ad96c0f04b60aa9a782c4de11356c2e06 Mon Sep 17 00:00:00 2001 From: helenarichie Date: Sat, 13 Jan 2024 17:24:46 -0500 Subject: [PATCH] get floor values from paramter file and remove density floor kernel --- builds/make.type.dust | 3 +++ src/global/global.cpp | 12 +++++++++++ src/global/global.h | 13 ++++++++---- src/grid/grid3D.cpp | 18 ++++++++++------ src/grid/grid3D.h | 3 ++- src/grid/initial_conditions.cpp | 15 +++++++------ src/hydro/hydro_cuda.cu | 35 ++----------------------------- src/hydro/hydro_cuda.h | 2 -- src/hydro/hydro_cuda_tests.cu | 1 - src/integrators/VL_3D_cuda.cu | 9 ++------ src/integrators/VL_3D_cuda.h | 2 +- src/integrators/simple_3D_cuda.cu | 4 ++-- src/integrators/simple_3D_cuda.h | 2 +- 13 files changed, 53 insertions(+), 66 deletions(-) diff --git a/builds/make.type.dust b/builds/make.type.dust index a6cf1bdd0..24a765302 100644 --- a/builds/make.type.dust +++ b/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 diff --git a/src/global/global.cpp b/src/global/global.cpp index 89347cbda..d0f898739 100644 --- a/src/global/global.cpp +++ b/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); diff --git a/src/global/global.h b/src/global/global.h index f7030b563..6c9524001 100644 --- 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 @@ -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 diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp index a76417321..d844784d5 100644 --- a/src/grid/grid3D.cpp +++ b/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 { diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index aff94c898..5354acf0a 100644 --- 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 index 87c5d392f..0e5aabc90 100644 --- a/src/grid/initial_conditions.cpp +++ b/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,17 +1339,17 @@ 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; @@ -1357,8 +1357,7 @@ void Grid3D::Clouds() 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 } } } diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu index e50cdd3f0..125e851b4 100644 --- a/src/hydro/hydro_cuda.cu +++ b/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; } } diff --git a/src/hydro/hydro_cuda.h b/src/hydro/hydro_cuda.h index 92a997790..8fcfbba05 100644 --- a/src/hydro/hydro_cuda.h +++ b/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); diff --git a/src/hydro/hydro_cuda_tests.cu b/src/hydro/hydro_cuda_tests.cu index 0403478e5..0796a3064 100644 --- a/src/hydro/hydro_cuda_tests.cu +++ b/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); diff --git a/src/integrators/VL_3D_cuda.cu b/src/integrators/VL_3D_cuda.cu index f299eda3c..25462fa91 100644 --- a/src/integrators/VL_3D_cuda.cu +++ b/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 diff --git a/src/integrators/VL_3D_cuda.h b/src/integrators/VL_3D_cuda.h index 3f2cf8d75..4104493bc 100644 --- a/src/integrators/VL_3D_cuda.h +++ b/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(); diff --git a/src/integrators/simple_3D_cuda.cu b/src/integrators/simple_3D_cuda.cu index c572dcd97..d3beb14b2 100644 --- a/src/integrators/simple_3D_cuda.cu +++ b/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 diff --git a/src/integrators/simple_3D_cuda.h b/src/integrators/simple_3D_cuda.h index 585c553ba..e2cea247e 100644 --- a/src/integrators/simple_3D_cuda.h +++ b/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();