From e2affb4398116f670408f64a990f26fcdc20aaba Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Wed, 13 Sep 2023 14:44:16 -0400 Subject: [PATCH 1/2] Add MHD support to slice outputs --- src/io/io.cpp | 120 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 111 insertions(+), 9 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 32fb8804d..c9817cc6d 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -15,6 +15,7 @@ #endif // HDF5 #include "../grid/grid3D.h" #include "../io/io.h" +#include "../utils/cuda_utilities.h" #include "../utils/timing_functions.h" // provides ScopedTimer #ifdef MPI_CHOLLA #include "../mpi/mpi_routines.h" @@ -1865,6 +1866,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real)); dataset_buffer_mz = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real)); dataset_buffer_E = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real)); + #ifdef MHD + std::vector dataset_buffer_magnetic_x(H.nx_real * H.ny_real); + std::vector dataset_buffer_magnetic_y(H.nx_real * H.ny_real); + std::vector dataset_buffer_magnetic_z(H.nx_real * H.ny_real); + #endif // MHD #ifdef DE dataset_buffer_GE = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real)); #endif @@ -1875,19 +1881,38 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) // Copy the xy slices to the memory buffers for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + zslice * H.nx * H.ny; + id = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice, H.nx, H.ny); buf_id = j + i * H.ny_real; + #ifdef MHD + int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, j + H.n_ghost, zslice, H.nx, H.ny); + int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost - 1, zslice, H.nx, H.ny); + int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice - 1, H.nx, H.ny); + #endif // MHD #ifdef MPI_CHOLLA // When there are multiple processes, check whether this slice is in // your domain if (zslice >= nz_local_start && zslice < nz_local_start + nz_local) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (zslice - nz_local_start + H.n_ghost) * H.nx * H.ny; - #endif // MPI_CHOLLA + id = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice - nz_local_start + H.n_ghost, H.nx, + H.ny); + #ifdef MHD + int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, j + H.n_ghost, + zslice - nz_local_start + H.n_ghost, H.nx, H.ny); + int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost - 1, + zslice - nz_local_start + H.n_ghost, H.nx, H.ny); + int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, + zslice - nz_local_start + H.n_ghost - 1, H.nx, H.ny); + #endif // MHD + #endif // MPI_CHOLLA dataset_buffer_d[buf_id] = C.density[id]; dataset_buffer_mx[buf_id] = C.momentum_x[id]; dataset_buffer_my[buf_id] = C.momentum_y[id]; dataset_buffer_mz[buf_id] = C.momentum_z[id]; dataset_buffer_E[buf_id] = C.Energy[id]; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]); + dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]); + dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]); + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = C.GasEnergy[id]; #endif @@ -1905,6 +1930,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my[buf_id] = 0; dataset_buffer_mz[buf_id] = 0; dataset_buffer_E[buf_id] = 0; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0; + dataset_buffer_magnetic_y[buf_id] = 0; + dataset_buffer_magnetic_z[buf_id] = 0; + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = 0; #endif @@ -1924,6 +1954,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_xy"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_xy"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_xy"); + #ifdef MHD + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_xy"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_xy"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_xy"); + #endif // MHD #ifdef DE status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_xy"); #endif @@ -1957,6 +1992,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real)); dataset_buffer_mz = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real)); dataset_buffer_E = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real)); + #ifdef MHD + dataset_buffer_magnetic_x.resize(H.nx_real * H.nz_real); + dataset_buffer_magnetic_y.resize(H.nx_real * H.nz_real); + dataset_buffer_magnetic_z.resize(H.nx_real * H.nz_real); + #endif // MHD #ifdef DE dataset_buffer_GE = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real)); #endif @@ -1967,19 +2007,38 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) // Copy the xz slices to the memory buffers for (k = 0; k < H.nz_real; k++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + yslice * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice, k + H.n_ghost, H.nx, H.ny); buf_id = k + i * H.nz_real; + #ifdef MHD + int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, yslice, k + H.n_ghost, H.nx, H.ny); + int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - 1, k + H.n_ghost, H.nx, H.ny); + int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice, k + H.n_ghost - 1, H.nx, H.ny); + #endif // MHD #ifdef MPI_CHOLLA // When there are multiple processes, check whether this slice is in // your domain if (yslice >= ny_local_start && yslice < ny_local_start + ny_local) { - id = (i + H.n_ghost) + (yslice - ny_local_start + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - #endif // MPI_CHOLLA + id = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost, k + H.n_ghost, H.nx, + H.ny); + #ifdef MHD + int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, yslice - ny_local_start + H.n_ghost, + k + H.n_ghost, H.nx, H.ny); + int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost - 1, + k + H.n_ghost, H.nx, H.ny); + int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost, + k + H.n_ghost - 1, H.nx, H.ny); + #endif // MHD + #endif // MPI_CHOLLA dataset_buffer_d[buf_id] = C.density[id]; dataset_buffer_mx[buf_id] = C.momentum_x[id]; dataset_buffer_my[buf_id] = C.momentum_y[id]; dataset_buffer_mz[buf_id] = C.momentum_z[id]; dataset_buffer_E[buf_id] = C.Energy[id]; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]); + dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]); + dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]); + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = C.GasEnergy[id]; #endif @@ -1997,6 +2056,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my[buf_id] = 0; dataset_buffer_mz[buf_id] = 0; dataset_buffer_E[buf_id] = 0; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0; + dataset_buffer_magnetic_y[buf_id] = 0; + dataset_buffer_magnetic_z[buf_id] = 0; + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = 0; #endif @@ -2016,6 +2080,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_xz"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_xz"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_xz"); + #ifdef MHD + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_xz"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_xz"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_xz"); + #endif // MHD #ifdef DE status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_xz"); #endif @@ -2050,6 +2119,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real)); dataset_buffer_mz = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real)); dataset_buffer_E = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real)); + #ifdef MHD + dataset_buffer_magnetic_x.resize(H.ny_real * H.nz_real); + dataset_buffer_magnetic_y.resize(H.ny_real * H.nz_real); + dataset_buffer_magnetic_z.resize(H.ny_real * H.nz_real); + #endif // MHD #ifdef DE dataset_buffer_GE = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real)); #endif @@ -2060,19 +2134,37 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) // Copy the yz slices to the memory buffers for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - id = xslice + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny); buf_id = k + j * H.nz_real; + #ifdef MHD + int id_xm1 = cuda_utilities::compute1DIndex(xslice - 1, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny); + int id_ym1 = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost - 1, k + H.n_ghost, H.nx, H.ny); + int id_zm1 = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost, k + H.n_ghost - 1, H.nx, H.ny); + #endif // MHD #ifdef MPI_CHOLLA // When there are multiple processes, check whether this slice is in // your domain if (xslice >= nx_local_start && xslice < nx_local_start + nx_local) { - id = (xslice - nx_local_start) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - #endif // MPI_CHOLLA + id = cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny); + #ifdef MHD + int id_xm1 = + cuda_utilities::compute1DIndex(xslice - nx_local_start - 1, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny); + int id_ym1 = + cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost - 1, k + H.n_ghost, H.nx, H.ny); + int id_zm1 = + cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost, k + H.n_ghost - 1, H.nx, H.ny); + #endif // MHD + #endif // MPI_CHOLLA dataset_buffer_d[buf_id] = C.density[id]; dataset_buffer_mx[buf_id] = C.momentum_x[id]; dataset_buffer_my[buf_id] = C.momentum_y[id]; dataset_buffer_mz[buf_id] = C.momentum_z[id]; dataset_buffer_E[buf_id] = C.Energy[id]; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]); + dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]); + dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]); + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = C.GasEnergy[id]; #endif @@ -2090,6 +2182,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) dataset_buffer_my[buf_id] = 0; dataset_buffer_mz[buf_id] = 0; dataset_buffer_E[buf_id] = 0; + #ifdef MHD + dataset_buffer_magnetic_x[buf_id] = 0; + dataset_buffer_magnetic_y[buf_id] = 0; + dataset_buffer_magnetic_z[buf_id] = 0; + #endif // MHD #ifdef DE dataset_buffer_GE[buf_id] = 0; #endif @@ -2109,6 +2206,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id) status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_yz"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_yz"); status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_yz"); + #ifdef MHD + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_yz"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_yz"); + status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_yz"); + #endif // MHD #ifdef DE status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_yz"); #endif From 9ab72fe7ae04f7e0fe4413f442191d473c58d004 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Wed, 13 Sep 2023 15:47:50 -0400 Subject: [PATCH 2/2] Move default paremeters into declaration Some member variables of the `parameters` struct were set to default values in an unrelated function. Moved the default initialization to the declaration of the struct. --- src/global/global.cpp | 16 ---------------- src/global/global.h | 23 ++++++++++++----------- 2 files changed, 12 insertions(+), 27 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index a4c697d3c..2cdb508d1 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -131,22 +131,6 @@ void parse_params(char *param_file, struct parameters *parms, int argc, char **a exit(1); return; } - // set default hydro file output parameter - parms->n_hydro = 1; - parms->n_particle = 1; - parms->n_slice = 1; - parms->n_projection = 1; - parms->n_rotated_projection = 1; - -#ifdef ROTATED_PROJECTION - // initialize rotation parameters to zero - parms->delta = 0; - parms->theta = 0; - parms->phi = 0; - parms->n_delta = 0; - parms->ddelta_dt = 0; - parms->flag_delta = 0; -#endif /*ROTATED_PROJECTION*/ #ifdef COSMOLOGY // Initialize file name as an empty string diff --git a/src/global/global.h b/src/global/global.h index b037c931d..75fee01fb 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -184,11 +184,11 @@ struct parameters { Real gamma; char init[MAXLEN]; int nfile; - int n_hydro; - int n_particle; - int n_projection; - int n_rotated_projection; - int n_slice; + int n_hydro = 1; + int n_particle = 1; + int n_projection = 1; + int n_rotated_projection = 1; + int n_slice = 1; int n_out_float32 = 0; int out_float32_density = 0; int out_float32_momentum_x = 0; @@ -275,16 +275,17 @@ struct parameters { char snr_filename[MAXLEN]; #endif #ifdef ROTATED_PROJECTION + // initialize rotation parameters to zero int nxr; int nzr; - Real delta; - Real theta; - Real phi; + Real delta = 0; + Real theta = 0; + Real phi = 0; Real Lx; Real Lz; - int n_delta; - Real ddelta_dt; - int flag_delta; + int n_delta = 0; + Real ddelta_dt = 0; + int flag_delta = 0; #endif /*ROTATED_PROJECTION*/ #ifdef COSMOLOGY Real H0;