Commit 6f889ed8 authored by kar's avatar kar
Browse files

added support for non-regular grids to QSP

parent c904026c
......@@ -258,7 +258,9 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
const char T1_name, const char S1_name, const int NS,
const int NT, double T1_max, double S1_max, double T1_min,
double S1_min, const int Nx, const int Ny, const int Nz,
string filename, const int plotnum);
string filename, const int plotnum, bool grab_grids,
TArrayn::DTArray *xgrid, TArrayn::DTArray *ygrid,
TArrayn::DTArray *zgrid);
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
......
......@@ -26,7 +26,9 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
const char T1_name, const char S1_name, const int NS,
const int NT, double T1_max, double S1_max, double T1_min,
double S1_min, const int Nx, const int Ny, const int Nz,
string filename, const int plotnum) {
string filename, const int plotnum, bool grab_grids,
TArrayn::DTArray *xgrid, TArrayn::DTArray *ygrid,
TArrayn::DTArray *zgrid) {
int local_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
......@@ -175,7 +177,7 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
double hS_inv = 1 / hS;
double hT_inv = 1 / hT;
int *local_hist = (int *)calloc(NS * NT, sizeof(int));
double *local_hist = (double *)calloc(NS * NT, sizeof(double));
if (!local_hist) {
std::cout << "Bad memory allocation. Exiting QSPCount" << std::endl;
return;
......@@ -242,21 +244,37 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
} else if (idxS >= NS) {
idxS = 0;
}
local_hist[index(idxS, idxT, NS, NT)] += 1;
double volume_weight;
if (grab_grids) {
double xgrid_weight, ygrid_weight, zgrid_weight;
xgrid_weight = (*xgrid)(i, j, k);
if (Ny > 1) { // Only dereference if not NULL
ygrid_weight = (*ygrid)(i, j, k);
} else {
ygrid_weight = 1.0; // Be passive if not included
}
zgrid_weight = (*zgrid)(i, j, k);
volume_weight = xgrid_weight * ygrid_weight * zgrid_weight;
} else {
volume_weight = 1.0;
}
local_hist[index(idxS, idxT, NS, NT)] += volume_weight;
}
}
}
MPI_Barrier(MPI_COMM_WORLD); // Wait for everyone to finish
if (local_rank == 0) {
int *glob_hist = (int *)calloc(NS * NT, sizeof(int));
double *glob_hist = (double *)calloc(NS * NT, sizeof(double));
if (!glob_hist) {
std::cout << "Bad memory allocation. Exiting QSPCount" << std::endl;
free(local_hist);
return;
}
MPI_Reduce(local_hist, glob_hist, // send and receive buffers
NS * NT, MPI_INT, // count and datatype
NS * NT, MPI_DOUBLE, // count and datatype
MPI_SUM, 0, // Reduction operator and root process #
MPI_COMM_WORLD); // Communicator
filename = filename + "." + boost::lexical_cast<string>(plotnum) + ".csv";
......@@ -279,7 +297,7 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
free(glob_hist);
} else {
MPI_Reduce(local_hist, NULL, // send and receive buffers
NS * NT, MPI_INT, // count and datatype
NS * NT, MPI_DOUBLE, // count and datatype
MPI_SUM, 0, // Reduction operator and root process #
MPI_COMM_WORLD); // Communicator
}
......
......@@ -10,6 +10,7 @@
// Defines limits of intrinsic types. Used as default values for
// T1_max, T1_min, S1_max and/orS1_min
#include <cassert>
#include <limits>
// Tensor variables for indexing
......@@ -38,6 +39,7 @@ double T1_max, S1_max, T1_min, S1_min;
char T1_name, S1_name;
int NT, NS;
string QSP_filename;
bool grab_grids;
// physical parameters
double visco; // viscosity (m^2/s)
......@@ -76,6 +78,7 @@ class userControl : public BaseCase {
DTArray *temp1, *temp2, *temp3; // arrays for vortex stretching / enstrophy production
DTArray *temp4; // array for storing temporary array for baroclinic vorticity
DTArray *A11, *A12, *A13, *A22, *A23, *A33; //Arrays for components of S^2+Omega^2
DTArray *xgrid, *ygrid, *zgrid; // Arrays for storing grid data
/* Size of domain */
double length_x() const { return Lx; }
......@@ -129,10 +132,27 @@ class userControl : public BaseCase {
split(deriv_filenames.c_str(), ' ', fields); // populate that vector
if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or
do_Q_and_R or do_lambda2 or do_hist ) {
temp1 = alloc_array(Nx,Ny,Nz);
temp2 = alloc_array(Nx,Ny,Nz);
if ( do_enst_stretch )
temp3 = alloc_array(Nx,Ny,Nz);
// If user asks to grab the grids, we allocate arrays to store
// them in memory
if (grab_grids) {
xgrid = alloc_array(Nx, Ny, Nz);
if (Ny > 1) {
ygrid = alloc_array(Nx, Ny, Nz);
}
zgrid = alloc_array(Nx, Ny, Nz);
} else {
// If user doesn't want to grab grids, we make sure not to
// allocate arrays for them and to set the pointers to NULL.
xgrid = NULL;
ygrid = NULL;
zgrid = NULL;
}
if ( do_enst_stretch ) temp3 = alloc_array(Nx,Ny,Nz);
}
if ( do_barvor ) {
......@@ -414,13 +434,27 @@ class userControl : public BaseCase {
S1_name == 'T') {
init_tracer_restart("t", *temp1);
}
QSPCount(*temp1, u, v, w, T1_name, S1_name, NS, NT,
T1_max, S1_max, T1_min, S1_min,
Nx, Ny, Nz, QSP_filename, plotnum);
// If user asked to grab the grids, we populate the grids
// with the correct data from disk
if (grab_grids) {
do_mapping(*xgrid, *ygrid, *zgrid);
} else {
// Make sure that if the user didn't want us to grab the
// grids then we haven't allocated data to store them!
assert(xgrid == NULL);
assert(ygrid == NULL);
assert(zgrid == NULL);
}
QSPCount(*temp1, u, v, w, T1_name, S1_name, NS, NT, T1_max,
S1_max, T1_min, S1_min, Nx, Ny, Nz, QSP_filename,
plotnum, grab_grids, xgrid, ygrid, zgrid);
if (master()) {
fprintf(stdout, "Completed the write for QSP.%d\n", plotnum);
}
}
}
}
......@@ -493,6 +527,7 @@ int main(int argc, char ** argv) {
add_option("do_R",&do_R,false,"Calculate R?");
add_option("do_Q_and_R",&do_Q_and_R,false,"Calculate Q and R?");
add_option("do_lambda2",&do_lambda2,false,"Calculate Lambda2?");
add_option("grab_grids",&grab_grids,false,"Grab grid data? (NOTE: grid data is expected to be named xgrid, ygrid zgrid).");
add_option("do_hist",&do_hist,false,"Create QSP Data?");
add_option("T1",&T1_name,'t', "Name of tracer 1 for QSP. Valid values are t (for rho),u,v,w,T (for temp) or k for K.E.");
add_option("S1",&S1_name,'w', "Name of tracer 2 for QSP. Valid values are t (for rho),u,v,w,T (for temp) or k for K.E.");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment