Skip to content

Commit

Permalink
resolve merge conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
helenarichie committed Jan 14, 2024
2 parents d783bb7 + 5fc5560 commit 8923164
Show file tree
Hide file tree
Showing 14 changed files with 452 additions and 24 deletions.
2 changes: 1 addition & 1 deletion Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ DIRS := src src/analysis src/chemistry_gpu src/cooling src/cooling_grackle s
src/cpu src/global src/gravity src/gravity/paris src/grid src/hydro \
src/integrators src/io src/main.cpp src/main_tests.cpp src/mhd\
src/model src/mpi src/old_cholla src/particles src/reconstruction \
src/riemann_solvers src/system_tests src/utils src/dust
src/riemann_solvers src/system_tests src/utils src/dust src/cloud_tracking

SUFFIX ?= .$(TYPE).$(MACHINE)

Expand Down
4 changes: 3 additions & 1 deletion builds/make.type.dust
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ MPI_GPU ?=
DFLAGS += -DCUDA
DFLAGS += -DMPI_CHOLLA
DFLAGS += -DPRECISION=2
DFLAGS += -DPLMP
DFLAGS += -DPPMP
DFLAGS += -DHLLC

DFLAGS += -DDE
Expand All @@ -36,6 +36,8 @@ DFLAGS += -DCLOUDY_COOLING
DFLAGS += -DSLICES
DFLAGS += -DPROJECTION

DFLAGS += -DCLOUD_TRACKING

DFLAGS += $(OUTPUT)

DFLAGS += -DOUTPUT_ALWAYS
Expand Down
6 changes: 4 additions & 2 deletions builds/make.type.hydro
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,20 @@ DFLAGS += -DSIMPLE
#DFLAGS += -DVL

# Apply a density and temperature floor
DFLAGS += -DDENSITY_FLOOR
# DFLAGS += -DDENSITY_FLOOR
DFLAGS += -DTEMPERATURE_FLOOR

# Solve the Gas Internal Energy usisng a Dual Energy Formalism
#DFLAGS += -DDE

# Apply cooling on the GPU from precomputed tables
#DFLAGS += -DCOOLING_GPU
DFLAGS += -DCOOLING_GPU

# Measure the Timing of the different stages
#DFLAGS += -DCPU_TIME

#DFLAGS += -DCLOUD_TRACKING

# Select output format
# Can also add -DSLICES and -DPROJECTIONS
OUTPUT ?= -DOUTPUT -DHDF5
Expand Down
159 changes: 159 additions & 0 deletions src/cloud_tracking/cloud_tracking.cu
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
/*!
* \file cloud_tracking.cu
* \author Helena Richie ([email protected])
* \brief
*/

#ifdef CLOUD_TRACKING

// STL includes
#include <stdio.h>

#include <cstdio>
#include <fstream>
#include <vector>

// Local includes

#include "../cloud_tracking/cloud_tracking.h"
#include "../global/global.h"
#include "../global/global_cuda.h"
#include "../grid/grid3D.h"
#include "../grid/grid_enum.h"
#include "../utils/DeviceVector.h"
#include "../utils/cuda_utilities.h"
#include "../utils/gpu.hpp"
#include "../utils/hydro_utilities.h"
#include "../utils/reduction_utilities.h"

void Cloud_Frame_Update(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz, int n_ghost,
int n_fields, Real dt, Real gamma, Real density_cloud_init, Real *integrand,
Real *density_cloud_tot, Real *mass_cloud_tot)
{
// cuda_utilities::AutomaticLaunchParams static const launchParams(Cloud_Tracking_Kernel)

int n_cells = nx * ny * nz;
int ngrid = (n_cells + TPB - 1) / TPB;
dim3 dim1dGrid(ngrid, 1, 1);
dim3 dim1dBlock(TPB, 1, 1);

// Allocate the device memory
cuda_utilities::DeviceVector<Real> static integrand_cloud(1, true);
cuda_utilities::DeviceVector<Real> static density_cloud(1, true);
cuda_utilities::DeviceVector<Real> static mass_cloud(1, true);

cuda_utilities::AutomaticLaunchParams static const launchParams(Cloud_Tracking_Kernel);

// printf("TPB: %d\n", TPB);
// printf("ngrid: %d\n", ngrid);
// printf("nx, ny, nz: %d %d %d\n", nx, ny, nz);
// printf("ngrid: %d\n", n_cells);
// printf("dim1dGrid: %d\n", dim1dGrid);
// printf("dim1dBlock: %d\n", dim1dBlock);
// printf("numBLocks: %d\n", launchParams.numBlocks);
// printf("threadsPerBLock: %d\n", launchParams.threadsPerBlock);

hipLaunchKernelGGL(Cloud_Tracking_Kernel, launchParams.numBlocks, launchParams.threadsPerBlock, 0, 0, dev_conserved,
nx, ny, nz, dx, dy, dz, n_ghost, n_fields, dt, gamma, density_cloud_init, integrand_cloud.data(),
density_cloud.data(), mass_cloud.data());
GPU_Error_Check();

*integrand = integrand_cloud[0];
*density_cloud_tot = density_cloud[0];
*mass_cloud_tot = mass_cloud[0];
}

__global__ void Cloud_Tracking_Kernel(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz,
int n_ghost, int n_fields, Real dt, Real gamma, Real density_cloud_init,
Real *integrand_cloud, Real *density_cloud, Real *mass_cloud)
{
int xid, yid, zid, n_cells;
n_cells = nx * ny * nz;
__shared__ Real density_stride[TPB];
__shared__ Real velocity_x_stride[TPB];
__shared__ Real mass_stride[TPB];
Real density, velocity_x, mass;

for (int i = 0; i < TPB; i++) {
density_stride[i] = 0;
velocity_x_stride[i] = 0;
mass_stride[i] = 0;
}

// Grid stride loop to perform as much of the reduction as possible. The
// fact that `id` has type `size_t` is important. I'm not totally sure why
// but setting it to int results in some kind of silent over/underflow issue
// even though we're not hitting those kinds of numbers. Setting it to type
// uint or size_t fixes them
for (size_t id = threadIdx.x + blockIdx.x * blockDim.x; id < n_cells; id += blockDim.x * gridDim.x) {
// get a global thread ID
cuda_utilities::compute3DIndices(id, nx, ny, xid, yid, zid);

// 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 = dev_conserved[id + n_cells * grid_enum::density];
velocity_x = dev_conserved[id + n_cells * grid_enum::momentum_x] / density;
mass = density * dx * dy * dz;

if (density > (1 / 3 * (density_cloud_init / DENSITY_UNIT))) {
// printf("inside if statement: %d %e\n", id, mass);
density_stride[threadIdx.x] += density;
velocity_x_stride[threadIdx.x] += velocity_x;
mass_stride[threadIdx.x] += mass;
// Do grid-wide reduction to compute mass-averaged cloud velocity (Shin et al. 2008, eq. 9)
}
}
}
__syncthreads();

reduction_utilities::Grid_Reduction_Add(density_stride[threadIdx.x] * velocity_x_stride[threadIdx.x],
integrand_cloud);
reduction_utilities::Grid_Reduction_Add(density_stride[threadIdx.x], density_cloud);
reduction_utilities::Grid_Reduction_Add(mass_stride[threadIdx.x], mass_cloud);
}

void Update_Grid_Velocities(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
Real velocity_cloud, Real density_cloud_tot, Real mass_cloud_tot)
{
// cuda_utilities::AutomaticLaunchParams static const launchParams(Velocity_Update);
int n_cells = nx * ny * nz;
int ngrid = (n_cells + TPB - 1) / TPB;
dim3 dim1dGrid(ngrid, 1, 1);
dim3 dim1dBlock(TPB, 1, 1);

hipLaunchKernelGGL(Velocity_Update, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, dt,
gamma, velocity_cloud, density_cloud_tot, mass_cloud_tot);
GPU_Error_Check();;
}

__global__ void Velocity_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt,
Real gamma, Real velocity_cloud, Real density_cloud_tot, Real mass_cloud_tot)
{
// get grid indices
int n_cells = nx * ny * nz;
int is, ie, js, je, ks, ke;
cuda_utilities::Get_Real_Indices(n_ghost, nx, ny, nz, is, ie, js, je, ks, ke);
// get a global thread ID
int blockId = blockIdx.x + blockIdx.y * gridDim.x;
int id = threadIdx.x + blockId * blockDim.x;
int id_z = id / (nx * ny);
int id_y = (id - id_z * nx * ny) / nx;
int id_x = id - id_z * nx * ny - id_y * nx;

Real density, momentum_x, velocity_x, energy;

Real energy_kinetic_cloud = 0.5 * mass_cloud_tot * pow(velocity_cloud, 2);

// threads corresponding to real cells do the calculation
if (id_x >= is && id_x < ie && id_y >= js && id_y < je && id_z >= ks && id_z < ke) {
density = dev_conserved[id + n_cells * grid_enum::density];
momentum_x = dev_conserved[id + n_cells * grid_enum::momentum_x];
velocity_x = density * momentum_x;
energy = dev_conserved[id + n_cells * grid_enum::Energy];
dev_conserved[id + n_cells * grid_enum::momentum_x] = (velocity_x - velocity_cloud) * density;
dev_conserved[id + n_cells * grid_enum::Energy] = energy - energy_kinetic_cloud;
}
}

#endif // CLOUD_TRACKING
31 changes: 31 additions & 0 deletions src/cloud_tracking/cloud_tracking.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
/*!
* \file cloud_tracking.h
* \author Helena Richie ([email protected])
* \brief Contains declaration for the kernel that does frame of reference tracking to track clouds.
*
*/
#ifdef CLOUD_TRACKING
#ifndef CLOUD_TRACKING_CUDA_H
#define CLOUD_TRACKING_CUDA_H

#include <math.h>

#include "../global/global.h"
#include "../utils/gpu.hpp"

void Cloud_Frame_Update(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz, int n_ghost,
int n_fields, Real dt, Real gamma, Real density_cloud_init, Real *integrand,
Real *density_cloud_tot, Real *mass_cloud_tot);

__global__ void Cloud_Tracking_Kernel(Real *dev_conserved, int nx, int ny, int nz, Real dx, Real dy, Real dz,
int n_ghost, int n_fields, Real dt, Real gamma, Real density_cloud_init,
Real *integrand_cloud, Real *density_cloud, Real *mass_cloud);

void Update_Grid_Velocities(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
Real velocity_cloud, Real density_cloud_tot, Real mass_cloud_tot);

__global__ void Velocity_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt,
Real gamma, Real velocity_cloud, Real density_cloud_tot, Real mass_cloud_tot);

#endif // CLOUD_TRACKING_CUDA_H
#endif // CLOUD_TRACKING
13 changes: 9 additions & 4 deletions src/global/global.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -101,10 +101,11 @@ char *Trim(char *s)
}

// NOLINTNEXTLINE(cert-err58-cpp)
const std::set<const char *> optionalParams = {
"flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta",
"delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Init_redshift",
"End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"};
const std::set<const char *> optionalParams = {"flag_delta", "ddelta_dt", "n_delta", "Lz",
"Lx", "phi", "theta", "delta",
"nzr", "nxr", "H0", "Omega_M",
"Omega_L", "Init_redshift", "End_redshift", "tile_length",
"n_proc_x", "n_proc_y", "n_proc_z", "density_cloud_init"};

/*! \fn int Is_Param_Valid(char *name);
* \brief Verifies that a param is valid (even if not needed). Avoids
Expand Down Expand Up @@ -453,6 +454,10 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
} else if (strcmp(name, "skewersdir") == 0) {
strncpy(parms->skewersdir, value, MAXLEN);
#endif
#endif
#ifdef CLOUD_TRACKING
} else if (strcmp(name, "density_cloud_init") == 0) {
parms->density_cloud_init = atof(value);
#endif
} else if (!Is_Param_Valid(name)) {
chprintf("WARNING: %s/%s: Unknown parameter/value pair!\n", name, value);
Expand Down
4 changes: 3 additions & 1 deletion src/global/global.h
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,9 @@ struct Parameters {
#ifdef TILED_INITIAL_CONDITIONS
Real tile_length;
#endif // TILED_INITIAL_CONDITIONS

#ifdef CLOUD_TRACKING
Real density_cloud_init;
#endif
#ifdef SET_MPI_GRID
// Set the MPI Processes grid [n_proc_x, n_proc_y, n_proc_z]
int n_proc_x;
Expand Down
2 changes: 1 addition & 1 deletion src/grid/cuda_boundaries.cu
Original file line number Diff line number Diff line change
Expand Up @@ -318,7 +318,7 @@ __global__ void Wind_Boundary_kernel(Real *c_device, int nx, int ny, int nz, int
d_0 = n_0 * mu * MP / DENSITY_UNIT;
P_0 = n_0 * KB * T_0 / PRESSURE_UNIT;

vx = 100 * TIME_UNIT / KPC; // km/s * (cholla unit conversion)
vx = 1000 * TIME_UNIT / KPC; // km/s * (cholla unit conversion)
vy = 0.0;
vz = 0.0;

Expand Down
42 changes: 42 additions & 0 deletions src/grid/grid3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,9 @@
#ifdef CLOUDY_COOL
#include "../cooling/load_cloudy_texture.h" // provides Load_Cuda_Textures and Free_Cuda_Textures
#endif
#ifdef CLOUD_TRACKING
#include "../cloud_tracking/cloud_tracking.h"
#endif

#ifdef PARALLEL_OMP
#include "../utils/parallel_omp.h"
Expand Down Expand Up @@ -464,6 +467,26 @@ void Grid3D::Execute_Hydro_Integrator(void)
} else if (H.nx > 1 && H.ny > 1 && H.nz > 1) // 3D
{
#ifdef CUDA
#ifdef CLOUD_TRACKING
// ==Subtract average cloud velocity from grid==
// Variables to store the mass-averaged cloud velocity and total cloud mass
Real velocity_cloud;
Real integrand, density_cloud_tot, mass_cloud_tot;
Real density;
// Do the grid-wide reduction to get the sum of rho*vx and total density for the entire cloud
Cloud_Frame_Update(C.device, H.nx, H.ny, H.nz, H.dx, H.dy, H.dz, H.n_ghost, H.n_fields, H.dt, gama,
H.density_cloud_init, &integrand, &density_cloud_tot, &mass_cloud_tot);

// chprintf("Cloud mass = %e\n", mass_cloud_tot);
velocity_cloud = integrand / mass_cloud_tot;
Update_Grid_Velocities(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama, velocity_cloud,
density_cloud_tot, mass_cloud_tot);
// chprintf("Cloud frame update = %d\n", state);
chprintf("Average cloud velocity = %e\n", velocity_cloud);
chprintf("Integrand = %e\n", integrand);
chprintf("Mass cloud = %e\n", mass_cloud_tot);
#endif // CLOUD_TRACKING

#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,
Expand Down Expand Up @@ -527,6 +550,24 @@ Real Grid3D::Update_Hydro_Grid()
Dust_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama);
#endif // DUST

#ifdef CLOUD_TRACKING
// ==Subtract average cloud velocity from grid==
// Variables to store the mass-averaged cloud velocity and total cloud mass
Real velocity_cloud;
Real integrand, density_cloud_tot, mass_cloud_tot;
Real density;
// Do the grid-wide reduction to get the sum of rho*vx and total density for the entire cloud
Cloud_Frame_Update(C.device, H.nx, H.ny, H.nz, H.dx, H.dy, H.dz, H.n_ghost, H.n_fields, H.dt, gama,
H.density_cloud_init, &integrand, &density_cloud_tot, &mass_cloud_tot);
velocity_cloud = integrand / mass_cloud_tot;
Update_Grid_Velocities(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama, velocity_cloud,
density_cloud_tot, mass_cloud_tot);
// chprintf("Cloud frame update = %d\n", state);
chprintf("Average cloud velocity = %e\n", velocity_cloud);
chprintf("Integrand = %e\n", integrand);
chprintf("Mass cloud = %e\n", mass_cloud_tot);
#endif // CLOUD_TRACKING

#endif // CUDA

#ifdef CHEMISTRY_GPU
Expand All @@ -536,6 +577,7 @@ Real Grid3D::Update_Hydro_Grid()
Timer.Chemistry.RecordTime(Chem.H.runtime_chemistry_step);
non_hydro_elapsed_time += Chem.H.runtime_chemistry_step;
#endif

C.HI_density = &C.host[H.n_cells * grid_enum::HI_density];
C.HII_density = &C.host[H.n_cells * grid_enum::HII_density];
C.HeI_density = &C.host[H.n_cells * grid_enum::HeI_density];
Expand Down
6 changes: 5 additions & 1 deletion src/grid/grid3D.h
Original file line number Diff line number Diff line change
Expand Up @@ -214,6 +214,10 @@ struct Header {
Real min_dt_slow;
#endif

#ifdef CLOUD_TRACKING
Real density_cloud_init;
#endif // CLOUD_TRACKING

/*! \var t_wall
* \brief Wall time */
Real t_wall;
Expand Down Expand Up @@ -675,7 +679,7 @@ class Grid3D
* gravitational collapse */
void Spherical_Overdensity_3D();

void Clouds();
void Clouds(struct Parameters P);

void Uniform_Grid();

Expand Down
Loading

0 comments on commit 8923164

Please sign in to comment.