Skip to content
Open
Show file tree
Hide file tree
Changes from 8 commits
Commits
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
22 changes: 0 additions & 22 deletions src/global/global.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
13 changes: 0 additions & 13 deletions src/global/global.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
35 changes: 0 additions & 35 deletions src/grid/grid3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
74 changes: 7 additions & 67 deletions src/grid/grid3D.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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. */
Expand Down
143 changes: 143 additions & 0 deletions src/io/RotatedProjWriter.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
/*!
* \file
* Implements the RotatedProjWriter type
*/

#include "RotatedProjWriter.h"

#include <cmath> // M_PI (note: not guaranteed by the C++ standard)
#ifdef HDF5
#include <hdf5.h>
#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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do think it's a good idea to check that the values of these parameters are valid, at least for the basic things of being positive/nonzero (either in this PR or a future one).

Also, I think it's okay for the n_delta/ddelta_dt parameters to go unused if they're defined with the wrong flag_delta value. But if you're concerned about it being confused, I think printing a warning could also be a good option.

// x-dir pixels in projection
this->nx = pmap.value<int>("nxr");
// z-dir pixels in projection
this->nz = pmap.value<int>("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<double>("Lx"));
// z-dir physical size of projection
this->Lz = Real(pmap.value<double>("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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that a switch statement would be clearer if you wanted to change it (but I don't think it's necessary for this PR if you don't feel like it).

if (this->rot_info_.flag_delta == 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++;
}

} else if (this->rot_info_.flag_delta == 2) {
// case 2 -- outputting 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);
} 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, 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);

#else
printf("Output_Rotated_Projected_Data only defined for HDF5 writes.\n");
#endif
}
Loading