diff --git a/src/Science.hpp b/src/Science.hpp
index 89d4c884a0dc956dabc65d9630cf030f65cd980a..cd5c38f754ead297aca9c7e09f6a338d2d164096 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -253,14 +253,38 @@ void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
/* Creates and outputs data for plotting a bivariate histogram with NS*NT bins
* between tracers S1 and T1 to "filename.csv" (don't include file extension).
*/
-void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
- const TArrayn::DTArray &v, const TArrayn::DTArray &w,
- 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, bool mapped,
- TArrayn::DTArray *xgrid, TArrayn::DTArray *ygrid,
- TArrayn::DTArray *zgrid);
+struct QSPOptions {
+ int NS;
+ int NT;
+ string filename;
+ double S1_max;
+ double S1_min;
+ double T1_max;
+ double T1_min;
+ string T1_name;
+ string S1_name;
+};
+
+struct QSPData {
+ TArrayn::DTArray *u;
+ TArrayn::DTArray *v;
+ TArrayn::DTArray *w;
+ TArrayn::DTArray *temp;
+ TArrayn::DTArray *rho;
+ TArrayn::DTArray *salinity;
+ TArrayn::DTArray *custom_T1;
+ TArrayn::DTArray *custom_S1;
+ TArrayn::DTArray *xgrid;
+ TArrayn::DTArray *ygrid;
+ TArrayn::DTArray *zgrid;
+ int Nx;
+ int Ny;
+ int Nz;
+ int plotnum;
+ bool mapped;
+};
+
+void QSPCount(QSPOptions qsp_options, QSPData qsp_data);
// 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 e2f73e92bb97613f1a849e3312a7bc6e926bcded..25fdae08abca743c8d0c48e222fd7722e134f65a 100644
--- a/src/Science/QSPCount.cpp
+++ b/src/Science/QSPCount.cpp
@@ -58,175 +58,282 @@ public:
~QSPVector() { free(data); }
};
-int index(int row, int col, int rows, int cols) {
- assert(row >= 0 && row < rows);
- assert(col >= 0 && col < cols);
- return (row * cols) + col;
+enum QSPType {
+ QSP_u,
+ QSP_v,
+ QSP_w,
+ QSP_ke,
+ QSP_temp,
+ QSP_rho,
+ QSP_salinity,
+ QSP_custom,
+};
+
+void QSP_write(int local_rank, const QSPVector &local_hist,
+ const QSPOptions &qsp_options, const QSPData &qsp_data) {
+ if (local_rank == 0) {
+ QSPVector glob_hist(qsp_options.NS, qsp_options.NT);
+ MPI_Reduce(local_hist.raw(), glob_hist.raw(), // send and receive buffers
+ qsp_options.NS * qsp_options.NT, // Count
+ MPI_DOUBLE, // datatype
+ MPI_SUM, 0, // Reduction operator and root process #
+ MPI_COMM_WORLD); // Communicator
+ string filename = qsp_options.filename + "." +
+ boost::lexical_cast<string>(qsp_data.plotnum) + ".csv";
+ std::fstream outfile;
+ outfile.open(filename.c_str(), std::ios_base::out);
+ if (outfile.is_open()) {
+ outfile << qsp_options.T1_max << ',' << qsp_options.T1_min << ','
+ << qsp_options.S1_max << ',' << qsp_options.S1_min;
+ for (int i = 4; i < qsp_options.NT; i++) {
+ outfile << ',' << 0;
+ }
+ outfile << std::endl;
+ for (int ii = 0; ii < qsp_options.NS; ii++) {
+ outfile << glob_hist(ii, 0);
+ for (int jj = 1; jj < qsp_options.NT; jj++) {
+ outfile << ',' << glob_hist(ii, jj);
+ }
+ outfile << std::endl;
+ }
+ }
+ } else {
+ MPI_Reduce(local_hist.raw(), NULL, // send and receive buffers
+ qsp_options.NS * qsp_options.NT, // count
+ MPI_DOUBLE, // datatype
+ MPI_SUM, 0, // Reduction operator and root process
+ MPI_COMM_WORLD); // Communicator
+ }
}
-void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
- const TArrayn::DTArray &v, const TArrayn::DTArray &w,
- 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, bool mapped,
- TArrayn::DTArray *xgrid, TArrayn::DTArray *ygrid,
- TArrayn::DTArray *zgrid) {
+QSPType QSPConvert(const std::string &name) {
+ QSPType converted_type;
+ if (name.compare("u") == 0) {
+ converted_type = QSP_u;
+ } else if (name.compare("v") == 0) {
+ converted_type = QSP_v;
+ } else if (name.compare("w") == 0) {
+ converted_type = QSP_w;
+ } else if (name.compare("ke") == 0) {
+ converted_type = QSP_ke;
+ } else if (name.compare("temp") == 0) {
+ converted_type = QSP_temp;
+ } else if (name.compare("rho") == 0) {
+ converted_type = QSP_rho;
+ } else if (name.compare("salinity") == 0) {
+ converted_type = QSP_salinity;
+ } else if (name.compare("custom") == 0) {
+ converted_type = QSP_custom;
+ }
+ return converted_type;
+}
- int local_rank;
- MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
- const TArrayn::DTArray *T1_ptr = NULL, *S1_ptr = NULL;
+TArrayn::DTArray *QSPPtr(const QSPData &qsp_data, const QSPType &type) {
+ TArrayn::DTArray *ptr = NULL;
- switch (T1_name) {
- case 'u':
- T1_ptr = &u;
+ switch (type) {
+ case QSP_u:
+ ptr = qsp_data.u;
break;
- case 'v':
- T1_ptr = &v;
+ case QSP_v:
+ ptr = qsp_data.v;
break;
- case 'w':
- T1_ptr = &w;
+ case QSP_w:
+ ptr = qsp_data.w;
break;
- case 'k':
+ case QSP_ke: // This is an odd case.
break;
- case 't':
- T1_ptr = &t;
+ case QSP_rho:
+ ptr = qsp_data.rho;
break;
- case 'T':
- T1_ptr = &t;
+ case QSP_temp:
+ ptr = qsp_data.temp;
break;
- default:
- return;
- }
-
- switch (S1_name) {
- case 'u':
- S1_ptr = &u;
- break;
- case 'v':
- S1_ptr = &v;
- break;
- case 'w':
- S1_ptr = &w;
- break;
- case 'k':
+ case QSP_salinity:
+ ptr = qsp_data.salinity;
break;
- case 't':
- S1_ptr = &t;
+ case QSP_custom: // Make sure to set the pointer yourself. Can't do it here.
break;
- case 'T':
- S1_ptr = &t;
- break;
- default:
- return;
}
- int i_low = u.lbound(blitz::firstDim);
- int j_low = u.lbound(blitz::secondDim);
- int k_low = u.lbound(blitz::thirdDim);
- int i_high = u.ubound(blitz::firstDim);
- int j_high = u.ubound(blitz::secondDim);
- int k_high = u.ubound(blitz::thirdDim);
+ return ptr;
+}
- // This block calculates the global min/max values, in case the user
- // didn't want to specify them.
+// compute the minimum and maximum values for either tracer, then substitute
+// any default values of +- infinity with the global min/max values instead.
+void QSPMaxMin(const QSPType &T1_type, const QSPType &S1_type,
+ TArrayn::DTArray *T1_ptr, TArrayn::DTArray *S1_ptr,
+ QSPOptions &qsp_options, const QSPData &qsp_data, int i_low,
+ int j_low, int k_low, int i_high, int j_high, int k_high) {
double double_max = std::numeric_limits<double>::max();
- double double_min = std::numeric_limits<double>::min();
- if (T1_max == double_max || S1_max == double_max || S1_max == double_min ||
- S1_min == double_min) {
- // If One of the variables is K.E or rho, we need to hand-roll the max/min
- // since, in general, the index of max(u) is the same as the index of max(v)
- if (T1_name == 'k' || T1_name == 't' || S1_name == 'k' || S1_name == 't') {
- double ke_max = -double_max, ke_min = double_max;
- double rho_max = -double_max, rho_min = double_max;
-
- // Main hand-rolled loop
- for (int i = i_low; i <= i_high; i++) {
- for (int j = j_low; j <= j_high; j++) {
- for (int k = k_low; k <= k_high; k++) {
- double ke_current = 0, tmp = 0;
- if (Nx > 1) {
- tmp = u(i, j, k);
- ke_current += tmp * tmp;
- }
- if (Ny > 1) {
- tmp = v(i, j, k);
- ke_current += tmp * tmp;
- }
- if (Nz > 1) {
- tmp = w(i, j, k);
+ // If One of the variables is K.E or rho, we need to hand-roll the max / min
+ // since, in general, the index of max(u) is the same as the index of max(v)
+ if (T1_type == QSP_ke || T1_type == QSP_rho || S1_type == QSP_ke ||
+ S1_type == QSP_rho) {
+ double ke_max = -double_max, ke_min = double_max;
+ double rho_max = -double_max, rho_min = double_max;
+ // Main hand-rolled loop
+ for (int i = i_low; i <= i_high; i++) {
+ for (int j = j_low; j <= j_high; j++) {
+ for (int k = k_low; k <= k_high; k++) {
+ double tmp;
+ if (T1_type == QSP_ke || S1_type == QSP_ke) {
+ double ke_current = 0;
+ tmp = (*qsp_data.u)(i, j, k);
+ ke_current += tmp * tmp;
+ if (qsp_data.Ny > 1) {
+ tmp = (*qsp_data.v)(i, j, k);
ke_current += tmp * tmp;
}
+ tmp = (*qsp_data.w)(i, j, k);
+ ke_current += tmp * tmp;
ke_current = 0.5 * ke_current;
- double rho_current = eqn_of_state_t(t(i, j, k));
ke_max = ke_current > ke_max ? ke_current : ke_max;
ke_min = ke_current < ke_min ? ke_current : ke_min;
+ }
+ if (T1_type == QSP_rho || S1_type == QSP_rho) {
+ double rho_current;
+ if (qsp_data.rho) {
+ rho_current = (*qsp_data.rho)(i, j, k);
+ } else if (qsp_data.temp && !qsp_data.salinity) {
+ rho_current = eqn_of_state_t((*qsp_data.temp)(i, j, k));
+ } else if (qsp_data.temp && qsp_data.salinity) {
+ rho_current = eqn_of_state((*qsp_data.temp)(i, j, k),
+ (*qsp_data.salinity)(i, j, k));
+ } else {
+ rho_current = 0;
+ }
rho_max = rho_current > rho_max ? rho_current : rho_max;
rho_min = rho_current < rho_min ? rho_current : rho_min;
}
}
}
+ }
+ double glob_ke_max, glob_ke_min, glob_rho_max, glob_rho_min;
+ MPI_Allreduce(&ke_max, &glob_ke_max, 1, MPI_DOUBLE, MPI_MAX,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(&ke_min, &glob_ke_min, 1, MPI_DOUBLE, MPI_MIN,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(&rho_max, &glob_rho_max, 1, MPI_DOUBLE, MPI_MAX,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(&rho_min, &glob_rho_min, 1, MPI_DOUBLE, MPI_MIN,
+ MPI_COMM_WORLD);
+ switch (T1_type) {
+ case QSP_ke:
+ qsp_options.T1_max =
+ qsp_options.T1_max == double_max ? glob_ke_max : qsp_options.T1_max;
+ qsp_options.T1_min =
+ qsp_options.T1_min == -double_max ? glob_ke_min : qsp_options.T1_min;
+ break;
+ case QSP_rho:
+ qsp_options.T1_max =
+ qsp_options.T1_max == double_max ? glob_rho_max : qsp_options.T1_max;
+ qsp_options.T1_min =
+ qsp_options.T1_min == -double_max ? glob_rho_min : qsp_options.T1_min;
+ break;
+ default:
+ qsp_options.T1_max = qsp_options.T1_max == double_max
+ ? psmax(max(*T1_ptr))
+ : qsp_options.T1_max;
+ qsp_options.T1_min = qsp_options.T1_min == -double_max
+ ? psmin(min(*T1_ptr))
+ : qsp_options.T1_min;
+ break;
+ }
+ switch (S1_type) {
+ case QSP_ke:
+ qsp_options.S1_max =
+ qsp_options.S1_max == double_max ? glob_ke_max : qsp_options.S1_max;
+ qsp_options.S1_min =
+ qsp_options.S1_min == -double_max ? glob_ke_min : qsp_options.S1_min;
+ break;
+ case QSP_rho:
+ qsp_options.S1_max =
+ qsp_options.S1_max == double_max ? glob_rho_max : qsp_options.S1_max;
+ qsp_options.S1_min =
+ qsp_options.S1_min == -double_max ? glob_rho_min : qsp_options.S1_min;
+ break;
+ default:
+ qsp_options.S1_max = qsp_options.S1_max == double_max
+ ? psmax(max(*S1_ptr))
+ : qsp_options.S1_max;
+ qsp_options.S1_min = qsp_options.S1_min == -double_max
+ ? psmin(min(*S1_ptr))
+ : qsp_options.S1_min;
+ break;
+ }
+ } else { // !(cond1 || cond2) == !cond1 && !cond2
+ qsp_options.S1_max = psmax(max(*S1_ptr));
+ qsp_options.S1_min = psmin(min(*S1_ptr));
+ qsp_options.T1_max = psmax(max(*T1_ptr));
+ qsp_options.T1_min = psmin(min(*T1_ptr));
+ }
+}
- double glob_ke_max, glob_ke_min, glob_rho_max, glob_rho_min;
- MPI_Allreduce(&ke_max, &glob_ke_max, 1, MPI_DOUBLE, MPI_MAX,
- MPI_COMM_WORLD);
- MPI_Allreduce(&ke_min, &glob_ke_min, 1, MPI_DOUBLE, MPI_MIN,
- MPI_COMM_WORLD);
- MPI_Allreduce(&rho_max, &glob_rho_max, 1, MPI_DOUBLE, MPI_MAX,
- MPI_COMM_WORLD);
- MPI_Allreduce(&rho_min, &glob_rho_min, 1, MPI_DOUBLE, MPI_MIN,
- MPI_COMM_WORLD);
+void QSPCount(QSPOptions qsp_options, QSPData qsp_data) {
- switch (T1_name) {
- case 'k':
- T1_max = glob_ke_max;
- T1_min = glob_ke_min;
- break;
- case 't':
- T1_max = glob_rho_max;
- T1_min = glob_rho_min;
- break;
- default:
- T1_max = psmax(max(*T1_ptr));
- T1_min = psmin(min(*T1_ptr));
- break;
- }
+ int local_rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
- switch (S1_name) {
- case 'k':
- S1_max = glob_ke_max;
- S1_min = glob_ke_min;
- break;
- case 't':
- S1_max = glob_rho_max;
- S1_min = glob_rho_min;
- break;
- default:
- S1_max = psmax(max(*S1_ptr));
- S1_min = psmin(min(*S1_ptr));
- break;
- }
- } else { // !(cond1 || cond2) == !cond1 && !cond2
- S1_max = psmax(max(*S1_ptr));
- S1_min = psmin(min(*S1_ptr));
- T1_max = psmax(max(*T1_ptr));
- T1_min = psmin(min(*T1_ptr));
+ // Find out what
+ QSPType S1_type = QSPConvert(qsp_options.S1_name);
+ QSPType T1_type = QSPConvert(qsp_options.T1_name);
+ TArrayn::DTArray *S1_ptr = QSPPtr(qsp_data, S1_type);
+ TArrayn::DTArray *T1_ptr = QSPPtr(qsp_data, T1_type);
+
+ if (T1_type == QSP_custom) {
+ T1_ptr = qsp_data.custom_T1;
+ }
+
+ if (S1_type == QSP_custom) {
+ S1_ptr = qsp_data.custom_S1;
+ }
+
+ if ((!S1_ptr && S1_type != QSP_ke) || (!T1_ptr && T1_type != QSP_ke)) {
+ std::cout << "Not enough data was provided for the requested tracer. "
+ "Aborting...\n";
+ return;
+ }
+
+ int i_low, j_low, k_low, i_high, j_high, k_high;
+
+ {
+ TArrayn::DTArray *temp_ptr;
+ if (S1_ptr) { // If S1 is not ke
+ temp_ptr = S1_ptr;
+ } else { // If S1 is ke we know u must exist
+ temp_ptr = qsp_data.u;
}
+ i_low = temp_ptr->lbound(blitz::firstDim);
+ j_low = temp_ptr->lbound(blitz::secondDim);
+ k_low = temp_ptr->lbound(blitz::thirdDim);
+ i_high = temp_ptr->ubound(blitz::firstDim);
+ j_high = temp_ptr->ubound(blitz::secondDim);
+ k_high = temp_ptr->ubound(blitz::thirdDim);
+ }
+
+ double double_max = std::numeric_limits<double>::max();
+ if (qsp_options.T1_max == double_max || qsp_options.S1_max == double_max ||
+ qsp_options.T1_min == -double_max || qsp_options.S1_min == -double_max) {
+ QSPMaxMin(T1_type, S1_type, T1_ptr, S1_ptr, qsp_options, qsp_data, i_low,
+ j_low, k_low, i_high, j_high, k_high);
}
- double hS = (S1_max - S1_min) / (double)NS;
- double hT = (T1_max - T1_min) / (double)NT;
+ double hS =
+ (qsp_options.S1_max - qsp_options.S1_min) / (double)qsp_options.NS;
+ double hT =
+ (qsp_options.T1_max - qsp_options.T1_min) / (double)qsp_options.NT;
double hS_inv = 1 / hS;
double hT_inv = 1 / hT;
- QSPVector local_hist(NS, NT);
- QSPVector global_z_max(Nx, Ny);
- QSPVector global_z_min(Nx, Ny);
+ QSPVector local_hist(qsp_options.NS, qsp_options.NT);
+ QSPVector global_z_max(qsp_data.Nx, qsp_data.Ny);
+ QSPVector global_z_min(qsp_data.Nx, qsp_data.Ny);
// Find the range of Lz values per 2D-slice
- if (mapped) {
- QSPVector local_z_max(Nx, Ny);
- QSPVector local_z_min(Nx, Ny);
+ if (qsp_data.mapped) {
+ QSPVector local_z_max(qsp_data.Nx, qsp_data.Ny);
+ QSPVector local_z_min(qsp_data.Nx, qsp_data.Ny);
// We are slicing as if we are doing zgrid[i, j, :]
for (int ii = i_low; ii <= i_high; ii++) {
for (int jj = j_low; jj <= j_high; jj++) {
@@ -236,7 +343,7 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
double tmp_z_max = -tmp_z_min;
double tmp;
for (int kk = k_low; kk <= k_high; kk++) {
- tmp = (*zgrid)(ii, jj, kk);
+ tmp = (*qsp_data.zgrid)(ii, jj, kk);
if (tmp > tmp_z_max) {
tmp_z_max = tmp;
} else if (tmp < tmp_z_min) {
@@ -247,10 +354,12 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
local_z_min(ii, jj) = tmp_z_min;
}
}
- MPI_Allreduce(local_z_min.raw(), global_z_min.raw(), Nx * Ny, MPI_DOUBLE,
- MPI_MIN, MPI_COMM_WORLD);
- MPI_Allreduce(local_z_max.raw(), global_z_max.raw(), Nx * Ny, MPI_DOUBLE,
- MPI_MAX, MPI_COMM_WORLD);
+ MPI_Allreduce(local_z_min.raw(), global_z_min.raw(),
+ qsp_data.Nx * qsp_data.Ny, MPI_DOUBLE, MPI_MIN,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(local_z_max.raw(), global_z_max.raw(),
+ qsp_data.Nx * qsp_data.Ny, MPI_DOUBLE, MPI_MAX,
+ MPI_COMM_WORLD);
}
// Main loop for QSP
@@ -259,85 +368,78 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
for (int j = j_low; j <= j_high; j++) {
for (int k = k_low; k <= k_high; k++) {
- if (T1_name == 'k') {
+ switch (T1_type) {
+ case QSP_ke:
Tval = 0;
- if (Nx > 1) {
- tmp = u(i, j, k);
- Tval += tmp * tmp;
- }
- if (Ny > 1) {
- tmp = v(i, j, k);
- Tval += tmp * tmp;
- }
- if (Nz > 1) {
- tmp = w(i, j, k);
+ tmp = (*qsp_data.u)(i, j, k);
+ Tval += tmp * tmp;
+ tmp = (*qsp_data.w)(i, j, k);
+ Tval += tmp * tmp;
+ if (qsp_data.Ny > 1) {
+ tmp = (*qsp_data.v)(i, j, k);
Tval += tmp * tmp;
}
Tval = 0.5 * Tval;
- } else if (T1_name == 't') {
+ break;
+ case QSP_rho:
tmp = (*T1_ptr)(i, j, k);
Tval = eqn_of_state_t(tmp);
- } else {
+ break;
+ default:
Tval = (*T1_ptr)(i, j, k);
- }
- int idxT = floor((Tval - T1_min) * hT_inv);
- if (idxT < 0) {
- idxT = 0;
- } else if (idxT >= NT) {
- idxT = 0;
+ break;
}
- if (S1_name == 'k') {
+ switch (S1_type) {
+ case QSP_ke:
Sval = 0;
- if (Nx > 1) {
- tmp = u(i, j, k);
- Sval += tmp * tmp;
- }
- if (Ny > 1) {
- tmp = v(i, j, k);
- Sval += tmp * tmp;
- }
- if (Nz > 1) {
- tmp = w(i, j, k);
+ tmp = (*qsp_data.u)(i, j, k);
+ Sval += tmp * tmp;
+ tmp = (*qsp_data.w)(i, j, k);
+ Sval += tmp * tmp;
+ if (qsp_data.Ny > 1) {
+ tmp = (*qsp_data.v)(i, j, k);
Sval += tmp * tmp;
}
Sval = 0.5 * Sval;
- } else if (S1_name == 't') {
+ break;
+ case QSP_rho:
tmp = (*S1_ptr)(i, j, k);
Sval = eqn_of_state_t(tmp);
- } else {
+ break;
+ default:
Sval = (*S1_ptr)(i, j, k);
+ break;
}
- int idxS = floor((Sval - S1_min) * hS_inv);
- if (idxS < 0) {
- idxS = 0;
- } else if (idxS >= NS) {
- idxS = 0;
- }
+
+ int idxS = floor((Sval - qsp_options.S1_min) * hS_inv);
+ int idxT = floor((Tval - qsp_options.T1_min) * hT_inv);
+ idxS = std::max(std::min(idxS, qsp_options.NS - 1), 0);
+ idxT = std::max(std::min(idxT, qsp_options.NT - 1), 0);
double volume_weight;
- if (mapped) {
+ if (qsp_data.mapped) {
// Calculate the Lz range
double Lzmax_now = global_z_max(i, j);
double Lzmin_now = global_z_min(i, j);
// Calculate the arc length
double arc, z_high, z_low, cos_high, cos_low;
- if (k > 0 && k < Nz - 1) {
+ if (k > 0 && k < qsp_data.Nz - 1) {
z_high = (double)(k + 1);
z_low = (double)(k - 1);
} else if (k == 0) {
z_high = (double)(k + 1);
z_low = (double)(k);
- } else if (k == Nz - 1) {
+ } else if (k == qsp_data.Nz - 1) {
z_high = (double)(k);
z_low = (double)(k - 1);
} else { // Failure
std::cout << "k was out of bounds somehow. Failing...\n";
return;
}
- cos_high = std::cos(M_PI * z_high / (double)Nz);
- cos_low = std::cos(M_PI * z_low / (double)Nz);
+ cos_high = std::cos(M_PI * z_high / (double)qsp_data.Nz);
+ cos_low = std::cos(M_PI * z_low / (double)qsp_data.Nz);
arc = 0.5 * (cos_low - cos_high);
// Calculate the volume weight
@@ -346,40 +448,11 @@ void QSPCount(const TArrayn::DTArray &t, const TArrayn::DTArray &u,
volume_weight = 1.0;
}
- // local_hist[index(idxS, idxT, NS, NT)] += volume_weight;
local_hist(idxS, idxT) += volume_weight;
}
}
}
MPI_Barrier(MPI_COMM_WORLD); // Wait for everyone to finish
- if (local_rank == 0) {
- QSPVector glob_hist(NS, NT);
- MPI_Reduce(local_hist.raw(), glob_hist.raw(), // send and receive buffers
- 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";
- std::fstream outfile;
- outfile.open(filename.c_str(), std::ios_base::out);
- if (outfile.is_open()) {
- outfile << T1_max << ',' << T1_min << ',' << S1_max << ',' << S1_min;
- for (int i = 4; i < NT; i++) {
- outfile << ',' << 0;
- }
- outfile << std::endl;
- for (int ii = 0; ii < NS; ii++) {
- outfile << glob_hist[index(ii, 0, NS, NT)];
- for (int jj = 1; jj < NT; jj++) {
- outfile << ',' << glob_hist[index(ii, jj, NS, NT)];
- }
- outfile << std::endl;
- }
- }
- } else {
- MPI_Reduce(local_hist.raw(), NULL, // send and receive buffers
- NS * NT, MPI_DOUBLE, // count and datatype
- MPI_SUM, 0, // Reduction operator and root process #
- MPI_COMM_WORLD); // Communicator
- }
+ QSP_write(local_rank, local_hist, qsp_options, qsp_data);
}
diff --git a/src/cases/derivatives/derivatives.cpp b/src/cases/derivatives/derivatives.cpp
index e4a962af4207f110fa56fa18f339170dd57869c3..9fd3b0db1004198ee3f364d79674f5f6e57c2254 100644
--- a/src/cases/derivatives/derivatives.cpp
+++ b/src/cases/derivatives/derivatives.cpp
@@ -36,7 +36,7 @@ const int z_ind = 2;
// QSP variables
double T1_max, S1_max, T1_min, S1_min;
-char T1_name, S1_name;
+string T1_name, S1_name;
int NT, NS;
string QSP_filename;
bool use_salinity;
@@ -152,7 +152,7 @@ class userControl : public BaseCase {
zgrid = NULL;
}
- if ( do_enst_stretch ) temp3 = alloc_array(Nx,Ny,Nz);
+ if ( do_enst_stretch or do_hist ) temp3 = alloc_array(Nx,Ny,Nz);
}
if ( do_barvor ) {
@@ -430,13 +430,6 @@ class userControl : public BaseCase {
// nor to make deep copies of them
if ( do_hist ){
- if (use_salinity) {
- init_tracer_restart("s", *temp1);
- } else if (T1_name == 't' || S1_name == 't' || T1_name == 'T' ||
- S1_name == 'T') {
- init_tracer_restart("t", *temp1);
- }
-
// If user asked to grab the grids, we populate the grids
// with the correct data from disk
if (mapped) {
@@ -449,9 +442,49 @@ class userControl : public BaseCase {
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, mapped, xgrid, ygrid, zgrid);
+ QSPOptions qsp_opts;
+ qsp_opts.S1_name = S1_name;
+ qsp_opts.T1_name = T1_name;
+ qsp_opts.filename = QSP_filename;
+ qsp_opts.NS = NS;
+ qsp_opts.NT = NT;
+ qsp_opts.S1_max = S1_max;
+ qsp_opts.S1_min = S1_min;
+ qsp_opts.T1_max = T1_max;
+ qsp_opts.T1_min = T1_min;
+
+ QSPData qsp_data;
+ qsp_data.mapped = mapped;
+ qsp_data.Nx = Nx;
+ qsp_data.Ny = Ny;
+ qsp_data.Nz = Nz;
+ qsp_data.plotnum = plotnum;
+ qsp_data.u = &u;
+ qsp_data.v = &v;
+ qsp_data.w = &w;
+ qsp_data.xgrid = xgrid;
+ qsp_data.ygrid = ygrid;
+ qsp_data.zgrid = zgrid;
+
+ qsp_data.temp = NULL;
+ qsp_data.rho = NULL;
+ qsp_data.salinity = NULL;
+
+ if (use_salinity) {
+ init_tracer_restart("s", *temp1);
+ qsp_data.salinity = temp1;
+ }
+ if (T1_name.compare("temp") == 0 || S1_name.compare("temp") == 0) {
+ init_tracer_restart("t", *temp2);
+ qsp_data.temp = temp2;
+ }
+ if (T1_name.compare("rho") == 0 || S1_name.compare("rho") == 0) {
+ init_tracer_restart("rho", *temp3);
+ qsp_data.rho = temp3;
+ }
+
+ QSPCount(qsp_opts, qsp_data);
+
if (master()) {
fprintf(stdout, "Completed the write for QSP.%d\n", plotnum);
}
@@ -530,13 +563,13 @@ int main(int argc, char ** argv) {
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("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.");
+ add_option("T1",&T1_name,"u", "Name of tracer 1 for QSP. Valid values are rho,u,v,w,temp or ke");
+ add_option("S1",&S1_name,"w", "Name of tracer 2 for QSP. Valid values are rho,u,v,w,temp or ke");
add_option("T1_max",&T1_max,std::numeric_limits<double>::max(), "Maximum explicit bin for T1 in QSP.");
add_option("T1_min",&T1_min,std::numeric_limits<double>::min(), "Minimum explicit bin for T1 in QSP.");
add_option("S1_max",&S1_max,std::numeric_limits<double>::max(), "Maximum explicit bin for S1 in QSP.");
add_option("S1_min",&S1_min,std::numeric_limits<double>::min(), "Minimum explicit bin for S1 in QSP.");
- add_option("QSP_salinity",&use_salinity, false, "Should salinity be read in from filename s?.");
+ add_option("salinity",&use_salinity, false, "Should salinity be read in from filename s?.");
add_option("QSP_filename",&QSP_filename,"QSP_default", "Filename to save data to. Don't include file extension.");
add_option("NS",&NS,10,"Number of bins for tracer S");
add_option("NT",&NT,10,"Number of bins for tracer T");
diff --git a/src/cases/derivatives/spins.conf b/src/cases/derivatives/spins.conf
index 9e9304c8f827a3a96740c154769aa2182d8c29a5..3d05b5ba8bed408714ff6da43f738ac245dee918 100644
--- a/src/cases/derivatives/spins.conf
+++ b/src/cases/derivatives/spins.conf
@@ -42,11 +42,11 @@ do_Q = false
do_R = false
do_Q_and_R = false
do_lambda2 = false
-do_hist = false
# QSP Parameters
-T1 = t
-S1 = k
+do_hist = false
+T1 = temp
+S1 = ke
QSP_filename = QSP
NT = 10
NS = 10
diff --git a/src/cases/qsp/qsp.cpp b/src/cases/qsp/qsp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3dca0407e408fa5b77c63b6296929c857f70efbc
--- /dev/null
+++ b/src/cases/qsp/qsp.cpp
@@ -0,0 +1,305 @@
+/* Derivative case file for computing derivatives of existing fields */
+// only able to return first and second derivatives
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // contains default class
+#include "../../Options.hpp" // config-file parser
+#include "../../Science.hpp" // Science content
+
+#include <cassert>
+#include <limits>
+#include <string>
+
+// Tensor variables for indexing
+blitz::firstIndex ii;
+blitz::secondIndex jj;
+blitz::thirdIndex kk;
+
+/* ------------------ Define parameters --------------------- */
+
+// Grid scales
+double Lx, Ly, Lz; // Grid lengths (m)
+int Nx, Ny, Nz; // Number of points in x, y, z
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+// Expansion types and indices
+S_EXP expan[3];
+const int x_ind = 0;
+const int y_ind = 1;
+const int z_ind = 2;
+
+// QSP variables
+double T1_max, S1_max, T1_min, S1_min;
+std::string T1_name, S1_name;
+int NT, NS;
+std::string QSP_filename;
+bool use_salinity, read_rho;
+std::string salinity_filename, temp_filename, rho_filename;
+std::string custom_T1_filename, custom_S1_filename;
+const double double_max = std::numeric_limits<double>::max();
+
+// current output number
+int plotnum;
+
+// Derivative options
+int start_sequence; // output number to start taking derivatives at
+int final_sequence; // output number to stop taking derivatives at
+int step_sequence; // step between outputs to take derivatives
+ // streching/tilting?
+bool v_exist; // Does the v field exist?
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+public:
+ /* Initialize things */
+ DTArray *xgrid, *ygrid, *zgrid; // Arrays for storing grid data
+
+ /* Size of domain */
+ double length_x() const { return Lx; }
+ double length_y() const { return Ly; }
+ double length_z() const { return Lz; }
+
+ /* Resolution in X, Y, and Z */
+ int size_x() const { return Nx; }
+ int size_y() const { return Ny; }
+ int size_z() const { return Nz; }
+
+ /* Set expansions (FREE_SLIP, NO_SLIP (in vertical) or PERIODIC) */
+ DIMTYPE type_x() const { return intype_x; }
+ DIMTYPE type_y() const { return intype_y; }
+ DIMTYPE type_z() const { return intype_z; }
+
+ /* Set other things */
+ int get_restart_sequence() const { return plotnum; }
+
+ /* Read grid (if mapped) */
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray &xg, DTArray &yg, DTArray &zg) {
+ init_grid_restart("x", "xgrid", xg);
+ if (Ny > 1)
+ init_grid_restart("y", "ygrid", yg);
+ init_grid_restart("z", "zgrid", zg);
+ }
+
+ /* Read fields and do derivatives */
+ void init_vels(DTArray &u, DTArray &v, DTArray &w) {
+ // If user asks to grab the grids, we allocate arrays to store
+ // them in memory
+ if (mapped) {
+ xgrid = alloc_array(Nx, Ny, Nz);
+ if (Ny > 1) {
+ ygrid = alloc_array(Nx, Ny, Nz);
+ }
+ zgrid = alloc_array(Nx, Ny, Nz);
+ do_mapping(*xgrid, *ygrid, *zgrid);
+ } 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;
+ }
+ QSPOptions qsp_opts;
+ qsp_opts.S1_name = S1_name;
+ qsp_opts.T1_name = T1_name;
+ qsp_opts.filename = QSP_filename;
+ qsp_opts.NS = NS;
+ qsp_opts.NT = NT;
+ qsp_opts.S1_max = S1_max;
+ qsp_opts.S1_min = S1_min;
+ qsp_opts.T1_max = T1_max;
+ qsp_opts.T1_min = T1_min;
+
+ QSPData qsp_data;
+ qsp_data.mapped = mapped;
+ qsp_data.Nx = Nx;
+ qsp_data.Ny = Ny;
+ qsp_data.Nz = Nz;
+ qsp_data.xgrid = xgrid;
+ qsp_data.ygrid = ygrid;
+ qsp_data.zgrid = zgrid;
+
+ // Compute derivatives at each requested output
+ for (plotnum = start_sequence; plotnum <= final_sequence;
+ plotnum = plotnum + step_sequence) {
+
+ // read in fields (if not already stored in memory)
+ // Compute QSP data. The code promises to not mutate the arrays,
+ // nor to make deep copies of them
+
+ // u
+ init_tracer_restart("u", u);
+ // v
+ if (v_exist) {
+ init_tracer_restart("v", v);
+ } else {
+ if (master())
+ fprintf(stdout, "No v field, setting v=0\n");
+ v = 0;
+ }
+ // w
+ init_tracer_restart("w", w);
+
+ qsp_data.plotnum = plotnum;
+ qsp_data.u = &u;
+ qsp_data.v = &v;
+ qsp_data.w = &w;
+ qsp_data.temp = NULL;
+ qsp_data.rho = NULL;
+ qsp_data.salinity = NULL;
+ qsp_data.custom_T1 = NULL;
+ qsp_data.custom_S1 = NULL;
+
+ DTArray *arr_salinity, *arr_temperature, *arr_rho, *custom_T1, *custom_S1;
+
+ if (use_salinity) {
+ arr_salinity = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(salinity_filename, *arr_salinity);
+ qsp_data.salinity = arr_salinity;
+ }
+
+ if (read_rho) {
+ arr_rho = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(rho_filename, *arr_rho);
+ qsp_data.rho = arr_rho;
+ }
+
+ if (T1_name.compare("temp") == 0 || S1_name.compare("temp") == 0) {
+ arr_temperature = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(temp_filename, *arr_temperature);
+ qsp_data.temp = arr_temperature;
+ }
+
+ if (T1_name.compare("custom") == 0) {
+ custom_T1 = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(custom_T1_filename, *custom_T1);
+ qsp_data.custom_T1 = custom_T1;
+ }
+
+ if (S1_name.compare("custom") == 0) {
+ custom_S1 = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(custom_S1_filename, *custom_S1);
+ qsp_data.custom_S1 = custom_S1;
+ }
+
+ QSPCount(qsp_opts, qsp_data);
+
+ if (master()) {
+ fprintf(stdout, "Completed QSP calculation for plotnum %d\n", plotnum);
+ }
+ }
+ }
+
+ // Constructor: Initialize local variables
+ userControl() {
+ compute_quadweights(size_x(), size_y(), size_z(), length_x(), length_y(),
+ length_z(), type_x(), type_y(), type_z());
+ }
+};
+
+/* The ''main'' routine */
+int main(int argc, char **argv) {
+ /* Initialize MPI. This is required even for single-processor runs,
+ since the inner routines assume some degree of parallelization,
+ even if it is trivial. */
+ MPI_Init(&argc, &argv);
+
+ /* ------------------ Define parameters from spins.conf ---------------------
+ */
+ options_init();
+
+ option_category("Grid Options");
+ add_option("Lx", &Lx, "Length of tank");
+ add_option("Ly", &Ly, 1.0, "Width of tank");
+ add_option("Lz", &Lz, "Height of tank");
+ add_option("Nx", &Nx, "Number of points in X");
+ add_option("Ny", &Ny, 1, "Number of points in Y");
+ add_option("Nz", &Nz, "Number of points in Z");
+
+ option_category("Grid mapping options");
+ add_option("mapped_grid", &mapped, false, "Is the grid mapped?");
+
+ string xgrid_type, ygrid_type, zgrid_type;
+ add_option("type_x", &xgrid_type,
+ "Grid type in X. Valid values are:\n"
+ " FOURIER: Periodic\n"
+ " FREE_SLIP: Cosine expansion\n"
+ " NO_SLIP: Chebyhsev expansion");
+ add_option("type_y", &ygrid_type, "FOURIER", "Grid type in Y");
+ add_option("type_z", &zgrid_type, "Grid type in Z");
+
+ option_category("QSP options");
+ add_option("start_sequence", &start_sequence,
+ "Sequence number to start taking derivatives at");
+ add_option("final_sequence", &final_sequence,
+ "Sequence number to stop taking derivatives at");
+ add_option("step_sequence", &step_sequence, 1,
+ "Step between outputs to take derivatives");
+ add_option("salinity_filename", &salinity_filename,
+ "Base Filename of Salinity data (optional).");
+ add_option("temp_filename", &temp_filename,
+ "Base Filename of temperature data (optional).");
+ add_option("T1_filename", &custom_T1_filename,
+ "Base Filename of custom tracer for T1.");
+ add_option("S1_filename", &custom_S1_filename,
+ "Base Filename of custom tracer for S1.");
+ add_option("T1", &T1_name, "u",
+ "Name of tracer 1 for QSP. Valid values are rho,u,v,w,temp, "
+ "salinity or ke");
+ add_option("S1", &S1_name, "w",
+ "Name of tracer 2 for QSP. Valid values are rho,u,v,w,temp, "
+ "salinity or ke");
+ add_option("T1_max", &T1_max, double_max,
+ "Maximum explicit bin for T1 in QSP.");
+ add_option("T1_min", &T1_min, -double_max,
+ "Minimum explicit bin for T1 in QSP.");
+ add_option("S1_max", &S1_max, double_max,
+ "Maximum explicit bin for S1 in QSP.");
+ add_option("S1_min", &S1_min, -double_max,
+ "Minimum explicit bin for S1 in QSP.");
+ add_option("salinity", &use_salinity, false,
+ "Should salinity be read in from a file?.");
+ add_option("read_rho", &read_rho, false,
+ "Should rho be read in from a file?.");
+ add_option("QSP_filename", &QSP_filename, "QSP_default",
+ "Filename to save data to. Don't include file extension.");
+ add_option("NS", &NS, 10, "Number of bins for tracer S");
+ add_option("NT", &NT, 10, "Number of bins for tracer T");
+ add_option("v_exist", &v_exist, "Does the v field exist?");
+ // Parse the options from the command line and config file
+ options_parse(argc, argv);
+
+ /* ------------------ Adjust and check parameters --------------------- */
+ /* Now, make sense of the options received. Many of these values
+ can be directly used, but the ones of string-type need further procesing.
+ */
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x,
+ intype_y, intype_z);
+ // vector of string types
+ grid_type[x_ind] = xgrid_type;
+ grid_type[y_ind] = ygrid_type;
+ grid_type[z_ind] = zgrid_type;
+
+ // adjust Ly for 2D
+ if (Ny == 1 and Ly != 1.0) {
+ Ly = 1.0;
+ if (master())
+ fprintf(stdout, "Simulation is 2 dimensional, "
+ "Ly has been changed to 1.0 for normalization.\n");
+ }
+
+ /* ------------------ Do stuff --------------------- */
+ userControl mycode; // Create an instantiated object of the above class
+ FluidEvolve<userControl> kevin_kh(&mycode);
+ kevin_kh.initialize();
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/qsp/spins.conf b/src/cases/qsp/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..dd65a3bd8b7742946d05410f553b805950966b6c
--- /dev/null
+++ b/src/cases/qsp/spins.conf
@@ -0,0 +1,69 @@
+## QSP configuration file
+
+# Spatial Parameters
+Lx = 6.4
+Ly = 1.0
+Lz = 0.3
+Nx = 4096
+Ny = 1
+Nz = 384
+
+# Expansion types
+type_x = FREE_SLIP
+type_y = FOURIER
+type_z = NO_SLIP
+mapped_grid = false
+v_exist = false
+
+# If salinity is needed (directly or indirectly) set this to true
+salinity = false
+salinity_filename = s
+
+# If you already pre-computed rho (the equation of state) and saved it to a
+# file, then set this to true. If this is false, and QSP can find the
+# temperature data but not rho, it will use the equation of state function from
+# Science.hpp, and if it also finds salinity it will add that to the equation
+# of state calculation too.
+read_rho = false
+rho_filename = rho
+
+# Sequence Options
+start_sequence = 0
+final_sequence = 0
+step_sequence = 1
+
+# Choose what tracers you want to use. T1 and S1 have no distinguishable
+# meaning, just stay consistent with NT and NS. Do note that if you pick ke
+# then QSP will compute it in-place. If you already pre-computed ke and have a
+# saved file for it, use the custom property explained below instead.
+# Tracer options: rho, temp, salinity, u, v, w, ke
+#
+# Note: If you need to use a custom tracer (i.e. not one of the above), set the
+# tracer name to "custom." If you do this however, it is your responsability to
+# ensure that:
+# a) you provide a valid custom filename too and
+# b) all prepocessing is already completed and the data in the file is exactly
+# what you want to use
+T1 = temp
+S1 = ke
+
+# Tracer base filenames. These are only relevent parameters IFF you set the
+# according tracer name to "custom" above. Otherwise, these are ignored.
+T1_filename = vorticity
+S1_filename = vorticity
+
+# Set the min/max values of the first/last bin in each dimension respectively.
+# These are optional parameters. If you remove these lines from the config QSP
+# will find the global min/max values for each timestep and use that instead.
+# So if you need fixed values for these make sure to include them!
+T1_min = 1
+T1_max = 1
+S1_min = 1
+S1_max = 1
+
+# Number of bins for T1 and S1 chosen above.
+NT = 10
+NS = 10
+
+# Base filename to save results
+QSP_filename = QSP