From 22fe5af47c65404ade9452c07267c7b94f725a6b Mon Sep 17 00:00:00 2001 From: Alwin Date: Mon, 14 Aug 2023 15:36:44 -0400 Subject: [PATCH] add velocity and temperature ceilings --- src/grid/grid3D.cpp | 15 +++++ src/hydro/hydro_cuda.cu | 120 ++++++++++++++++++++++++++++++++++++++++ src/hydro/hydro_cuda.h | 14 ++++- 3 files changed, 146 insertions(+), 3 deletions(-) diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp index 1777a4d81..5fd4ac7d8 100644 --- a/src/grid/grid3D.cpp +++ b/src/grid/grid3D.cpp @@ -499,6 +499,21 @@ Real Grid3D::Update_Grid(void) #endif #endif + #ifdef VELOCITY_CEILING + const Real V_ceiling_cholla = 0.005; // roughly 10000 km/s + Velocity_Ceiling(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, gama, V_ceiling_cholla); +#endif //VELOCITY_CEILING + + // Temperature Ceiling + #ifdef TEMPERATURE_CEILING + // 1e51 ergs / (m_p * (pc/cm)^3) = 45000 km/s + // sqrt(1e10 K * kB/ m_mp) = 9000 km/s + const Real T_ceiling_kelvin = 5e9;// 1e10; + Temperature_Ceiling(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, gama, T_ceiling_kelvin); +#endif //TEMPERATURE_CEILING + + + #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell Real max_dti_slow; diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu index eb50990de..da4c09552 100644 --- a/src/hydro/hydro_cuda.cu +++ b/src/hydro/hydro_cuda.cu @@ -594,6 +594,126 @@ Real Calc_dt_GPU(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n return dev_dti[0]; } + +#ifdef VELOCITY_CEILING + +__global__ void Velocity_Ceiling_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real V_ceiling, int* counter) +{ + const int id = threadIdx.x + blockIdx.x * blockDim.x; + const int n_cells = nx * ny * nz; + int xid, yid, zid; + cuda_utilities::compute3DIndices(id, nx, ny, xid, yid, zid); + const bool real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost && zid > n_ghost - 1 && zid < nz - n_ghost); + if (!real_cell) return; + + const Real d = dev_conserved[id]; + const Real max_momentum = d * V_ceiling; + + const Real d_inv = 1.0 / d; + /* + const Real vx = dev_conserved[1 * n_cells + id] * d_inv; + const Real vy = dev_conserved[2 * n_cells + id] * d_inv; + const Real vz = dev_conserved[3 * n_cells + id] * d_inv; + const Real E = dev_conserved[4 * n_cells + id]; + */ + for (int momentum_index = 1; momentum_index <= 3; momentum_index++) { + // Reduce momentum if velocity is too large + + const Real momentum = dev_conserved[momentum_index * n_cells + id]; + if (abs(momentum) > max_momentum) { + const Real new_momentum = max_momentum*momentum/abs(momentum); + const Real diff_energy = 0.5 * d_inv * ((momentum * momentum) - (max_momentum*max_momentum)); + // Write in the new momentum + dev_conserved[momentum_index * n_cells + id] = new_momentum; + // Thermalize the energy + #ifdef DE + dev_conserved[(n_fields - 1) * n_cells + id] += diff_energy; + #endif //DE + atomicAdd(counter, 1); + } + } + + +} + +void Velocity_Ceiling(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real V_ceiling) +{ + + int n_cells = nx * ny * nz; + int ngrid = (n_cells + TPB - 1) / TPB; + dim3 dim1dGrid(ngrid, 1, 1); + dim3 dim1dBlock(TPB, 1, 1); + + cuda_utilities::DeviceVector counter(1, true); + int* dev_counter = counter.data(); + + if (nx > 1 && ny > 1 && nz > 1) { // 3D + hipLaunchKernelGGL(Velocity_Ceiling_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, gamma, V_ceiling, dev_counter); + } + int host_counter = counter[0]; + if (host_counter > 0) { + printf("HYDRO WARNING: Velocity Ceiling applied to num_cells: %d \n", host_counter); + } +} + + +#endif + + +#ifdef TEMPERATURE_CEILING + +__global__ void Temperature_Ceiling_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real T_ceiling, int* counter) +{ + const int id = threadIdx.x + blockIdx.x * blockDim.x; + const int n_cells = nx * ny * nz; + int xid, yid, zid; + cuda_utilities::compute3DIndices(id, nx, ny, xid, yid, zid); + const bool real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost && zid > n_ghost - 1 && zid < nz - n_ghost); + if (!real_cell) return; + + const Real d = dev_conserved[id]; + const Real d_inv = 1.0 / d; + const Real vx = dev_conserved[1 * n_cells + id] * d_inv; + const Real vy = dev_conserved[2 * n_cells + id] * d_inv; + const Real vz = dev_conserved[3 * n_cells + id] * d_inv; + const Real E = dev_conserved[4 * n_cells + id]; + + const Real temperature_Kelvin = (gamma - 1) * (E - 0.5 * (vx*vx + vy*vy + vz*vz) * d) * ENERGY_UNIT / (d * DENSITY_UNIT / 0.6 / MP) / KB; + + if (temperature_Kelvin <= T_ceiling) return; + + const Real factor = T_ceiling / temperature_Kelvin; + + dev_conserved[4 * n_cells + id] *= factor; + #ifdef DE + dev_conserved[(n_fields - 1) * n_cells + id] *= factor; + #endif //DE + atomicAdd(counter, 1); +} + +void Temperature_Ceiling(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real T_ceiling) +{ + + int n_cells = nx * ny * nz; + int ngrid = (n_cells + TPB - 1) / TPB; + dim3 dim1dGrid(ngrid, 1, 1); + dim3 dim1dBlock(TPB, 1, 1); + + cuda_utilities::DeviceVector counter(1, true); + int* dev_counter = counter.data(); + + if (nx > 1 && ny > 1 && nz > 1) { // 3D + hipLaunchKernelGGL(Temperature_Ceiling_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, gamma, T_ceiling, dev_counter); + } + int host_counter = counter[0]; + if (host_counter > 0) { + printf("HYDRO WARNING: Temperature Ceiling applied to num_cells: %d \n", host_counter); + } +} + + +#endif //TEMPERATURE_CEILING + #ifdef AVERAGE_SLOW_CELLS void Average_Slow_Cells(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dx, Real dy, diff --git a/src/hydro/hydro_cuda.h b/src/hydro/hydro_cuda.h index 9b4403be0..c071aa33f 100644 --- a/src/hydro/hydro_cuda.h +++ b/src/hydro/hydro_cuda.h @@ -75,15 +75,23 @@ __global__ void Sync_Energies_2D(Real *dev_conserved, int nx, int ny, int n_ghos __global__ void Sync_Energies_3D(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real gamma, int n_fields); +#ifdef VELOCITY_CEILING +void Velocity_Ceiling(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real V_ceiling); +#endif // VELOCITY CEILING + +#ifdef TEMPERATURE_CEILING +void Temperature_Ceiling(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real gamma, Real T_ceiling); +#endif // TEMPERATURE CEILING + #ifdef AVERAGE_SLOW_CELLS void Average_Slow_Cells(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dx, Real dy, Real dz, Real gamma, Real max_dti_slow, Real xbound, Real ybound, Real zbound, int nx_offset, - int ny_offset, int nz_offset); + int ny_offset, int nz_offset); __global__ void Average_Slow_Cells_3D(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dx, - Real dy, Real dz, Real gamma, Real max_dti_slow, Real xbound, Real ybound, - Real zbound, int nx_offset, int ny_offset, int nz_offset); + Real dy, Real dz, Real gamma, Real max_dti_slow, Real xbound, Real ybound, Real zbound, + int nx_offset, int ny_offset, int nz_offset); #endif #ifdef TEMPERATURE_FLOOR