diff --git a/src/global/global.cpp b/src/global/global.cpp index 63094637c..9a30a3ae4 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -287,28 +287,6 @@ bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameter } else if (strcmp(name, "prng_seed") == 0) { parms->prng_seed = atoi(value); #endif // PARTICLES -#ifdef ROTATED_PROJECTION - } else if (strcmp(name, "nxr") == 0) { - parms->nxr = atoi(value); - } else if (strcmp(name, "nzr") == 0) { - parms->nzr = atoi(value); - } else if (strcmp(name, "delta") == 0) { - parms->delta = atof(value); - } else if (strcmp(name, "theta") == 0) { - parms->theta = atof(value); - } else if (strcmp(name, "phi") == 0) { - parms->phi = atof(value); - } else if (strcmp(name, "Lx") == 0) { - parms->Lx = atof(value); - } else if (strcmp(name, "Lz") == 0) { - parms->Lz = atof(value); - } else if (strcmp(name, "n_delta") == 0) { - parms->n_delta = atoi(value); - } else if (strcmp(name, "ddelta_dt") == 0) { - parms->ddelta_dt = atof(value); - } else if (strcmp(name, "flag_delta") == 0) { - parms->flag_delta = atoi(value); -#endif /*ROTATED_PROJECTION*/ #ifdef TILED_INITIAL_CONDITIONS } else if (strcmp(name, "tile_length") == 0) { parms->tile_length = atof(value); diff --git a/src/global/global.h b/src/global/global.h index ef9665c27..f6e98c556 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -292,19 +292,6 @@ struct Parameters { char sw_filename[MAXLEN]; #endif #endif -#ifdef ROTATED_PROJECTION - // initialize rotation parameters to zero - int nxr; - int nzr; - Real delta = 0; - Real theta = 0; - Real phi = 0; - Real Lx; - Real Lz; - int n_delta = 0; - Real ddelta_dt = 0; - int flag_delta = 0; -#endif /*ROTATED_PROJECTION*/ #ifdef COSMOLOGY Real H0; Real Omega_M; diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp index 34b30a4ed..6928360e5 100644 --- a/src/grid/grid3D.cpp +++ b/src/grid/grid3D.cpp @@ -213,41 +213,6 @@ void Grid3D::Initialize(struct Parameters *P) // allocate memory AllocateMemory(); -#ifdef ROTATED_PROJECTION - // x-dir pixels in projection - R.nx = P->nxr; - // z-dir pixels in projection - R.nz = P->nzr; - // minimum x location to project - R.nx_min = 0; - // minimum z location to project - R.nz_min = 0; - // maximum x location to project - R.nx_max = R.nx; - // maximum z location to project - R.nz_max = R.nz; - // rotation angle about z direction - R.delta = M_PI * (P->delta / 180.); // convert to radians - // rotation angle about x direction - R.theta = M_PI * (P->theta / 180.); // convert to radians - // rotation angle about y direction - R.phi = M_PI * (P->phi / 180.); // convert to radians - // x-dir physical size of projection - R.Lx = P->Lx; - // z-dir physical size of projection - R.Lz = P->Lz; - // initialize a counter for rotated outputs - R.i_delta = 0; - // number of rotated outputs in a complete revolution - R.n_delta = P->n_delta; - // rate of rotation between outputs, for an actual simulation - R.ddelta_dt = P->ddelta_dt; - // are we not rotating about z(0)? - // are we outputting multiple rotations(1)? or rotating during a - // simulation(2)? - R.flag_delta = P->flag_delta; -#endif /*ROTATED_PROJECTION*/ - // Values for lower limit for density and temperature #ifdef TEMPERATURE_FLOOR H.temperature_floor = P->temperature_floor; diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index 1cdaaab39..3e3aae3d6 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -50,67 +50,11 @@ #include "../analysis/analysis.h" #endif -struct Rotation { - /*! \var nx - * \brief Number of pixels in x-dir of rotated, projected image*/ - int nx; - - /*! \var nz - * \brief Number of pixels in z-dir of rotated, projected image*/ - int nz; - - /*! \var nx_min - * \brief Left most point in the projected image for this subvolume*/ - int nx_min; - - /*! \var nx_max - * \brief Right most point in the projected image for this subvolume*/ - int nx_max; - - /*! \var nz_min - * \brief Bottom most point in the projected image for this subvolume*/ - int nz_min; - - /*! \var nz_max - * \brief Top most point in the projected image for this subvolume*/ - int nz_max; - - /*! \var delta - * \brief Rotation angle about z axis in simulation frame*/ - Real delta; - - /*! \var theta - * \brief Rotation angle about x axis in simulation frame*/ - Real theta; - - /*! \var phi - * \brief Rotation angle about y axis in simulation frame*/ - Real phi; - - /*! \var Lx - * \brief Physical x-dir size of projected image*/ - Real Lx; - - /*! \var Lz - * \brief Physical z-dir size of projected image*/ - Real Lz; - - /*! \var i_delta - * \brief number of output projection for delta rotation*/ - int i_delta; - - /*! \var n_delta - * \brief total number of output projection for delta rotation*/ - Real n_delta; - - /*! \var ddelta_dt - * \brief rate of delta rotation*/ - Real ddelta_dt; - - /*! \var flag_delta - * \brief output mode for box rotation*/ - int flag_delta; -}; +// forward-declare the Rotation struct +namespace io +{ +struct Rotation; +} struct Header { /*! \var n_cells @@ -293,10 +237,6 @@ class Grid3D * \brief Header for the grid */ struct Header H; - /*! \var struct Rotation R - * \brief Rotation struct for data projections */ - struct Rotation R; - #ifdef GRAVITY // Object that contains data for gravity Grav3D Grav; @@ -508,12 +448,12 @@ class Grid3D /*! \fn void Write_Header_Rotated_HDF5(hid_t file_id) * \brief Write the relevant header info to the HDF5 file for rotated * projection. */ - void Write_Header_Rotated_HDF5(hid_t file_id); + void Write_Header_Rotated_HDF5(hid_t file_id, io::Rotation &R); /*! \fn void Write_Rotated_Projection_HDF5(hid_t file_id) * \brief Write rotated projected data to a file, at the current simulation * time. */ - void Write_Rotated_Projection_HDF5(hid_t file_id); + void Write_Rotated_Projection_HDF5(hid_t file_id, const io::Rotation &R); /*! \fn void Write_Slices_HDF5(hid_t file_id) * \brief Write xy, xz, and yz slices of all data to a file. */ diff --git a/src/io/RotatedProjWriter.cpp b/src/io/RotatedProjWriter.cpp new file mode 100644 index 000000000..83c8eb213 --- /dev/null +++ b/src/io/RotatedProjWriter.cpp @@ -0,0 +1,149 @@ +/*! + * \file + * Implements the RotatedProjWriter type + */ + +#include "RotatedProjWriter.h" + +#include // M_PI (note: not guaranteed by the C++ standard) +#ifdef HDF5 + #include +#endif // HDF5 + +#include "../global/global.h" // Parameters +#include "../grid/grid3D.h" +#include "../io/io.h" +#include "../utils/error_handling.h" + +io::Rotation::Rotation(ParameterMap &pmap) +{ + // thoughts for the future: + // do we want to enforce valid domains for arguments? We might require that + // - nxr & nzr are positive + // - delta, phi, & theta satisfy 0 <= x < 180 (maybe we want to be more flexible) + // - Lx & Lz are positive (or at least non-zero) + // - I think flag_delta should only be 0, 1, or 2 + // + // should we make it an error to specify the: + // - ddelta_dt parameter when flag_delta != 1 + // - n_delta parameter when flag_delta != 2 + + // x-dir pixels in projection + this->nx = pmap.value("nxr"); + // z-dir pixels in projection + this->nz = pmap.value("nzr"); + // minimum x location to project + this->nx_min = 0; + // minimum z location to project + this->nz_min = 0; + // maximum x location to project + this->nx_max = this->nx; + // maximum z location to project + this->nz_max = this->nz; + // rotation angle about z direction + this->delta = Real(M_PI * (pmap.value_or("delta", 0.0) / 180.)); // convert to radians + // rotation angle about x direction + this->theta = Real(M_PI * (pmap.value_or("theta", 0.0) / 180.)); // convert to radians + // rotation angle about y direction + this->phi = Real(M_PI * (pmap.value_or("phi", 0.0) / 180.)); // convert to radians + // x-dir physical size of projection + this->Lx = Real(pmap.value("Lx")); + // z-dir physical size of projection + this->Lz = Real(pmap.value("Lz")); + // initialize a counter for rotated outputs + this->i_delta = 0; + // number of rotated outputs in a complete revolution + this->n_delta = pmap.value_or("n_delta", 0); + CHOLLA_ASSERT(this->n_delta >= 0, "the \"n_delta\" parameter must not be negative"); + // rate of rotation between outputs, for an actual simulation + this->ddelta_dt = Real(pmap.value_or("ddelta_dt", 0.0)); + // are we not rotating about z(0)? + // are we outputting multiple rotations(1)? or rotating during a + // simulation(2)? + this->flag_delta = pmap.value_or("flag_delta", 0); +} + +void io::RotatedProjWriter::operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template) +{ +#ifdef HDF5 + hid_t file_id; + herr_t status; + + // create the filename + std::string filename = fname_template.format_fname(nfile, "_rot_proj"); + + // it may be a little more explicit to use a switch statement instead of if/elif/else + switch (this->rot_info_.flag_delta) { + case 1: { + // if flag_delta==1, then we are just outputting a + // bunch of rotations of the same snapshot + char fname[200]; + + for (int i_delta = 0; i_delta < this->rot_info_.n_delta; i_delta++) { + filename += "." + std::to_string(this->rot_info_.i_delta); + chprintf("Outputting rotated projection %s.\n", fname); + + // determine delta about z by output index + this->rot_info_.delta = 2.0 * M_PI * ((double)i_delta) / ((double)this->rot_info_.n_delta); + + // Create a new file + file_id = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + + // Write header (file attributes) + G.Write_Header_Rotated_HDF5(file_id, this->rot_info_); + + // Write the density and temperature projections to the output file + G.Write_Rotated_Projection_HDF5(file_id, this->rot_info_); + + // Close the file + status = H5Fclose(file_id); + + CHOLLA_ASSERT(status >= 0, "Rotated Projection: File write failed. ProcID: %d\n", procID); + + // iterate this->rot_info_.i_delta + this->rot_info_.i_delta++; + } + break; + } + case 2: { // outputing at a rotating delta + // rotation rate given in the parameter file + this->rot_info_.delta = fmod(nfile * this->rot_info_.ddelta_dt * 2.0 * M_PI, (2.0 * M_PI)); + + // Create a new file + file_id = H5Fcreate(filename.data(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + + // Write header (file attributes) + G.Write_Header_Rotated_HDF5(file_id, this->rot_info_); + + // Write the density and temperature projections to the output file + G.Write_Rotated_Projection_HDF5(file_id, this->rot_info_); + + // Close the file + status = H5Fclose(file_id); + break; + } + case 0: { // just output at the delta given in the parameter file + + // Create a new file + file_id = H5Fcreate(filename.data(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + + // Write header (file attributes) + G.Write_Header_Rotated_HDF5(file_id, this->rot_info_); + + // Write the density and temperature projections to the output file + G.Write_Rotated_Projection_HDF5(file_id, this->rot_info_); + + // Close the file + status = H5Fclose(file_id); + break; + } + default: + CHOLLA_ERROR("Invalid flag_delta: %d", this->rot_info_.flag_delta); + } + + CHOLLA_ASSERT(status >= 0, "Rotated Projection: File write failed. ProcID: %d\n", procID); + +#else + printf("Output_Rotated_Projected_Data only defined for HDF5 writes.\n"); +#endif +} \ No newline at end of file diff --git a/src/io/RotatedProjWriter.h b/src/io/RotatedProjWriter.h new file mode 100644 index 000000000..6a1c28be4 --- /dev/null +++ b/src/io/RotatedProjWriter.h @@ -0,0 +1,90 @@ +/*! + * \file + * Declares the RotatedProjWriter type + */ + +#pragma once + +#include "../global/global.h" +#include "../grid/grid3D.h" +#include "../io/FnameTemplate.h" // define FnameTemplate +#include "../io/ParameterMap.h" // define ParameterMap + +namespace io +{ + +/*! Tracks rotation information for creating Rotated Projections + */ +struct Rotation { + // the default constructor would put the instance in an invalid state + Rotation() = delete; + + /*! Construct a new instance */ + Rotation(ParameterMap &pmap); + + /*! Number of pixels in x-dir of rotated, projected image*/ + int nx; + + /*! Number of pixels in z-dir of rotated, projected image*/ + int nz; + + /*! Left most point in the projected image for this subvolume*/ + int nx_min; + + /*! Right most point in the projected image for this subvolume*/ + int nx_max; + + /*! Bottom most point in the projected image for this subvolume*/ + int nz_min; + + /*! Top most point in the projected image for this subvolume*/ + int nz_max; + + /*! Rotation angle about z axis in simulation frame*/ + Real delta; + + /*! Rotation angle about x axis in simulation frame*/ + Real theta; + + /*! Rotation angle about y axis in simulation frame*/ + Real phi; + + /*! Physical x-dir size of projected image*/ + Real Lx; + + /*! Physical z-dir size of projected image*/ + Real Lz; + + /*! number of output projection for delta rotation*/ + int i_delta; + + /*! total number of output projection for delta rotation*/ + Real n_delta; + + /*! rate of delta rotation*/ + Real ddelta_dt; + + /*! output mode for box rotation*/ + int flag_delta; +}; + +/*! \brief A callable that writes rotated projections + * + * \note + * The initial skeleton is basically a placeholder + */ +class RotatedProjWriter +{ + /*! Tracks the rotation information */ + Rotation rot_info_; + + public: + RotatedProjWriter() = delete; + RotatedProjWriter(ParameterMap &pmap) : rot_info_(pmap) {} + + /*! A callable method that writes a rotated projection of the grid data to file. + */ + void operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template); +}; + +} // namespace io \ No newline at end of file diff --git a/src/io/WriterManager.cpp b/src/io/WriterManager.cpp index 02631b317..3f5938352 100644 --- a/src/io/WriterManager.cpp +++ b/src/io/WriterManager.cpp @@ -10,7 +10,8 @@ #include #include "../gravity/grav3D.h" -#include "../io/ParameterMap.h" // define ParameterMap +#include "../io/ParameterMap.h" // ParameterMap +#include "../io/RotatedProjWriter.h" // RotatedProjWriter #include "../io/io.h" io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap) : fname_template_(P) @@ -36,8 +37,8 @@ io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap) : fnam #endif /*PROJECTION*/ #ifdef ROTATED_PROJECTION - packs_.push_back(io::detail::WriterPack{"rotated_projection", pmap.value_or("n_rotated_projection", 1), - &Output_Rotated_Projected_Data}); + packs_.push_back(io::detail::WriterPack{ + "rotated_projection", pmap.value_or("n_rotated_projection", 1), {io::RotatedProjWriter(pmap)}}); #endif /*ROTATED_PROJECTION*/ #ifdef SLICES diff --git a/src/io/io.cpp b/src/io/io.cpp index 93178fc09..6e9a21865 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -365,105 +365,6 @@ void Output_Projected_Data(Grid3D &G, struct Parameters P, int nfile, const Fnam #endif // HDF5 } -/* Output a rotated projection of the grid data to file. */ -void Output_Rotated_Projected_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template) -{ -#ifdef HDF5 - hid_t file_id; - herr_t status; - - // create the filename - std::string filename = fname_template.format_fname(nfile, "_rot_proj"); - - if (G.R.flag_delta == 1) { - // if flag_delta==1, then we are just outputting a - // bunch of rotations of the same snapshot - int i_delta; - char fname[200]; - - for (i_delta = 0; i_delta < G.R.n_delta; i_delta++) { - filename += "." + std::to_string(G.R.i_delta); - chprintf("Outputting rotated projection %s.\n", fname); - - // determine delta about z by output index - G.R.delta = 2.0 * M_PI * ((double)i_delta) / ((double)G.R.n_delta); - - // Create a new file - file_id = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - - // Write header (file attributes) - G.Write_Header_Rotated_HDF5(file_id); - - // Write the density and temperature projections to the output file - G.Write_Rotated_Projection_HDF5(file_id); - - // Close the file - status = H5Fclose(file_id); - #ifdef MPI_CHOLLA - if (status < 0) { - printf("Output_Rotated_Projected_Data: File write failed. ProcID: %d\n", procID); - chexit(-1); - } - #else - if (status < 0) { - printf("Output_Rotated_Projected_Data: File write failed.\n"); - exit(-1); - } - #endif - - // iterate G.R.i_delta - G.R.i_delta++; - } - - } else if (G.R.flag_delta == 2) { - // case 2 -- outputting at a rotating delta - // rotation rate given in the parameter file - G.R.delta = fmod(nfile * G.R.ddelta_dt * 2.0 * M_PI, (2.0 * M_PI)); - - // Create a new file - file_id = H5Fcreate(filename.data(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - - // Write header (file attributes) - G.Write_Header_Rotated_HDF5(file_id); - - // Write the density and temperature projections to the output file - G.Write_Rotated_Projection_HDF5(file_id); - - // Close the file - status = H5Fclose(file_id); - } else { - // case 0 -- just output at the delta given in the parameter file - - // Create a new file - file_id = H5Fcreate(filename.data(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - - // Write header (file attributes) - G.Write_Header_Rotated_HDF5(file_id); - - // Write the density and temperature projections to the output file - G.Write_Rotated_Projection_HDF5(file_id); - - // Close the file - status = H5Fclose(file_id); - } - - #ifdef MPI_CHOLLA - if (status < 0) { - printf("Output_Rotated_Projected_Data: File write failed. ProcID: %d\n", procID); - chexit(-1); - } - #else - if (status < 0) { - printf("Output_Rotated_Projected_Data: File write failed.\n"); - exit(-1); - } - #endif - -#else - printf("Output_Rotated_Projected_Data only defined for HDF5 writes.\n"); -#endif -} - /* Output xy, xz, and yz slices of the grid data. */ void Output_Slices(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template) @@ -689,7 +590,7 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) /*! \fn void Write_Header_Rotated_HDF5(hid_t file_id) * \brief Write the relevant header info to the HDF5 file for rotated * projection. */ -void Grid3D::Write_Header_Rotated_HDF5(hid_t file_id) +void Grid3D::Write_Header_Rotated_HDF5(hid_t file_id, io::Rotation &R) { hid_t attribute_id, dataspace_id; herr_t status; @@ -1622,7 +1523,7 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) /*! \fn void Write_Rotated_Projection_HDF5(hid_t file_id) * \brief Write rotated projected data to a file, at the current simulation * time. */ -void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) +void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id, const io::Rotation &R) { hid_t dataset_id, dataspace_xzr_id; Real *dataset_buffer_dxzr; diff --git a/src/io/io.h b/src/io/io.h index 32e03ff8b..360b312bb 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -7,6 +7,7 @@ #include "../global/global.h" #include "../grid/grid3D.h" #include "../io/FnameTemplate.h" +#include "../io/RotatedProjWriter.h" // io::Rotation #include "../io/WriterManager.h" /* Local function that designates whether we are using a root-process. It gives @@ -35,9 +36,6 @@ void Output_Float32(Grid3D& G, struct Parameters P, int nfile, const FnameTempla /* Output a projection of the grid data to file. */ void Output_Projected_Data(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template); -/* Output a rotated projection of the grid data to file. */ -void Output_Rotated_Projected_Data(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template); - /* Output xy, xz, and yz slices of the grid data to file. */ void Output_Slices(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template);