diff --git a/src/Science.hpp b/src/Science.hpp
index fb54033a5ee0319deec510c765a344621a100317..096d88f10d4abc8bd9f089ce59d7f6af4ea40dc8 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -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)
diff --git a/src/Science/QSPCount.cpp b/src/Science/QSPCount.cpp
index c2d2df524437e0fd2bea3d8485cfae37a9cfcc02..780763ae502b6145d8841cfc8b4bfd0e9c1cf0b0 100644
--- a/src/Science/QSPCount.cpp
+++ b/src/Science/QSPCount.cpp
@@ -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
}
diff --git a/src/cases/derivatives/derivatives.cpp b/src/cases/derivatives/derivatives.cpp
index 99b158abab73d374d9a00df52328ea5a1cbe9dcd..393d3b7a6b36e33030ac6fa29e114965e803431c 100644
--- a/src/cases/derivatives/derivatives.cpp
+++ b/src/cases/derivatives/derivatives.cpp
@@ -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.");