Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add MHD Support to Slice Outputs #334

Merged
merged 3 commits into from
Sep 14, 2023
Merged

Add MHD Support to Slice Outputs #334

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e2affb4
Add MHD support to slice outputs
bcaddy Sep 13, 2023
9ab72fe
Move default paremeters into declaration
bcaddy Sep 13, 2023
6609d9e
Merge branch 'dev' into dev-mhdSlicesOutput
evaneschneider Sep 14, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 0 additions & 16 deletions src/global/global.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
23 changes: 12 additions & 11 deletions src/global/global.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
120 changes: 111 additions & 9 deletions src/io/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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<Real> dataset_buffer_magnetic_x(H.nx_real * H.ny_real);
std::vector<Real> dataset_buffer_magnetic_y(H.nx_real * H.ny_real);
std::vector<Real> 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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down