diff --git a/matlab/spins_QSP_csv.m b/matlab/spins_QSP_csv.m
new file mode 100644
index 0000000000000000000000000000000000000000..63f60d07bd5011503a0587b9d266f0edfb3e7bd9
--- /dev/null
+++ b/matlab/spins_QSP_csv.m
@@ -0,0 +1,15 @@
+function [T1_max, T1_min, S1_max, S1_min, pdfCount] = spins_QSP_reader(filename)
+ % Reads in the csv file and interprets it in the correct way for the user.
+ %
+ % Input: Filename (string) name of csv file.
+ %
+ % Output: Returns the max/min values for T1 and S1 (as specified in
+ % spins.conf) and the (normalized) pdfCount histogram.
+ A = importdata(filename);
+ T1_max = A(1, 1);
+ T1_min = A(1, 2);
+ S1_max = A(1, 3);
+ S1_min = A(1, 4);
+ pdfCount = A(2:end, :);
+ pdfCount /= sum(sum(pdfCount));
+end
diff --git a/src/Science.hpp b/src/Science.hpp
index 24b459c0b039755d14cb37db5eec7fa9ead1b304..fb54033a5ee0319deec510c765a344621a100317 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -13,12 +13,12 @@
extern Array<double,1> _quadw_x, _quadw_y, _quadw_z;
// Marek's Overturning Diagnostic
-blitz::Array<double,3> overturning_2d(blitz::Array<double,3> const & rho,
+blitz::Array<double,3> overturning_2d(blitz::Array<double,3> const & rho,
blitz::Array<double,1> const & zgrid, TArrayn::Dimension reduce = TArrayn::thirdDim );
// Read in a 2D file and interpret it as a 2D slice of a 3D array, for
// initialization with read-in-data from a program like MATLAB
-void read_2d_slice(blitz::Array<double,3> & fillme, const char * filename,
+void read_2d_slice(blitz::Array<double,3> & fillme, const char * filename,
int Nx, int Ny);
void read_2d_restart(blitz::Array<double,3>& fillme, const char* filename,
@@ -40,7 +40,7 @@ void compute_baroclinic_vort_x(TArrayn::DTArray & barovortx, TArrayn::DTArray &
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_baroclinic_vort_y(TArrayn::DTArray & barovorty, TArrayn::DTArray & T,
TArrayn::Grad * gradient_op, const string * grid_type);
-void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp,
+void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp,
TArrayn::DTArray & T, TArrayn::Grad * gradient_op, const string * grid_type, bool v_exist);
// Enstrophy density
@@ -63,7 +63,7 @@ void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho = false);
// Quadrature weights
-void compute_quadweights(int szx, int szy, int szz,
+void compute_quadweights(int szx, int szy, int szz,
double Lx, double Ly, double Lz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
@@ -132,20 +132,20 @@ void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
inline double nleos_inline(double T, double S){
- // Returns the density (kg/m^3)
+ // Returns the density (kg/m^3)
// MacDougall et. al. 2003 (JAOT 20)
//
// This EOS is a rational function taken at pressure = 0
//
// Constants for Numerator
-
+
// First find any regions of negative salinity and set them to zero for this operation.
-
- if (S < 0) {
+
+ if (S < 0) {
S = 0;
};
- const double N0 = 9.99843699e+02;
+ const double N0 = 9.99843699e+02;
const double N1 = 7.35212840e+00;
const double N2 = -5.45928211e-02;
const double N3 = 3.98476704e-04;
@@ -162,7 +162,7 @@ inline double nleos_inline(double T, double S){
double numerator = N0 + T*(N1 + T*(N2 + N3*T)) + S*(N4 + N5*T + N6*S) + p*(N7 + N8*T*T + N9*S + p*(N10 + N11*T*T));
- // Constants for denominator
+ // Constants for denominator
const double D0 = 1.00000000e+00;
const double D1 = 7.28606739e-03;
const double D2 = -4.60835542e-05;
@@ -176,7 +176,7 @@ inline double nleos_inline(double T, double S){
const double D10 = 5.30848875e-06;
const double D11 = -3.03175128e-16;
const double D12 = -1.27934137e-17;
-
+
double denominator = D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(S)*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p));
return numerator/denominator;
@@ -187,12 +187,12 @@ void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
TArrayn::DTArray & S);
inline double compute_alpha(double T0, double S0){
- // Computes the thermal expansion coefficient at S0 and T0
+ // Computes the thermal expansion coefficient at S0 and T0
// Derivative is a finite difference for simplicity
// Does not divide by rho_0, that is done in lineos()
const double dT = 1e-08;
double TpdT = T0 + dT;
-
+
double alpha = (nleos_inline(TpdT,S0) - nleos_inline(T0,S0))/dT;
return alpha;
@@ -200,7 +200,7 @@ inline double compute_alpha(double T0, double S0){
BZ_DECLARE_FUNCTION(compute_alpha)
inline double compute_beta(double T0, double S0){
- // Computes the haline contraction coefficient at S0 and T0
+ // Computes the haline contraction coefficient at S0 and T0
// Derivative is a finite difference for simplicity
// Does not divide by rho_0, that is done in lineos()
const double dS = 1e-08;
@@ -211,9 +211,9 @@ inline double compute_beta(double T0, double S0){
return beta;
}
BZ_DECLARE_FUNCTION(compute_beta)
-
+
inline double compute_rho0(double T0, double S0){
- // Computes the reference density at S0 and T0
+ // Computes the reference density at S0 and T0
double rho_0 = nleos_inline(T0,S0);
return rho_0;
@@ -229,7 +229,7 @@ void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TA
/*
inline double fresh_quad(double T){
- // Returns the density (kg/m^3) for water using simple quadratic fit to
+ // Returns the density (kg/m^3) for water using simple quadratic fit to
// MacDougall et. al. 2003 (JAOT 20)
// Constants are determined via a quadratic fit preserving the Temperature of maximum density
const double rho_max = 9.999744074665388e+02; // Density of freshwater at Tmd (deg C)
@@ -248,7 +248,17 @@ void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type, TArrayn::DTArray & A11, TArrayn::DTArray & A12,
TArrayn::DTArray & A13, TArrayn::DTArray & A22, TArrayn::DTArray & A23,
- TArrayn::DTArray & A33);
+ TArrayn::DTArray & A33);
+
+/* 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);
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
@@ -283,7 +293,7 @@ inline double eqn_of_state_dT(double T){
// Constants are from table 4 of the above paper, at pressure 0
// (This is appropriate since this is currently an incompressible
// model)
-
+
// Kept the same names for these constants as in eqn_of_state
const double c2 = 5.10768e-2; // T term, constant after derivative
const double c4 = -7.40849e-3; // T^2 term, 2T after derivative
diff --git a/src/Science/QSPCount.cpp b/src/Science/QSPCount.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c2d2df524437e0fd2bea3d8485cfae37a9cfcc02
--- /dev/null
+++ b/src/Science/QSPCount.cpp
@@ -0,0 +1,287 @@
+#include "../Science.hpp"
+#include "boost/lexical_cast.hpp"
+#include <algorithm>
+#include <blitz/array.h>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <ios>
+#include <iostream>
+#include <limits>
+#include <mpi.h>
+#include <string>
+#include <vector>
+
+using namespace Transformer;
+using namespace blitz;
+
+int index(int row, int col, int rows, int cols) {
+ assert(row >= 0 && row < rows);
+ assert(col >= 0 && col < cols);
+ return (row * cols) + col;
+}
+
+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) {
+
+ int local_rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
+ const TArrayn::DTArray *T1_ptr = NULL, *S1_ptr = NULL;
+
+ switch (T1_name) {
+ case 'u':
+ T1_ptr = &u;
+ break;
+ case 'v':
+ T1_ptr = &v;
+ break;
+ case 'w':
+ T1_ptr = &w;
+ break;
+ case 'k':
+ break;
+ case 't':
+ T1_ptr = &t;
+ break;
+ case 'T':
+ T1_ptr = &t;
+ 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':
+ break;
+ case 't':
+ S1_ptr = &t;
+ 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);
+
+ // This block calculates the global min/max values, in case the user
+ // didn't want to specify them.
+ 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);
+ 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;
+ 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_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;
+ }
+
+ 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));
+ }
+ }
+
+ double hS = (S1_max - S1_min) / (double)NS;
+ double hT = (T1_max - T1_min) / (double)NT;
+ double hS_inv = 1 / hS;
+ double hT_inv = 1 / hT;
+
+ int *local_hist = (int *)calloc(NS * NT, sizeof(int));
+ if (!local_hist) {
+ std::cout << "Bad memory allocation. Exiting QSPCount" << std::endl;
+ return;
+ }
+
+ // Main loop for QSP
+ double Tval, Sval, tmp;
+ 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++) {
+
+ if (T1_name == 'k') {
+ 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);
+ Tval += tmp * tmp;
+ }
+ Tval = 0.5 * Tval;
+ } else if (T1_name == 't') {
+ tmp = (*T1_ptr)(i, j, k);
+ Tval = eqn_of_state_t(tmp);
+ } else {
+ 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;
+ }
+
+ if (S1_name == 'k') {
+ 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);
+ Sval += tmp * tmp;
+ }
+ Sval = 0.5 * Sval;
+ } else if (S1_name == 't') {
+ tmp = (*S1_ptr)(i, j, k);
+ Sval = eqn_of_state_t(tmp);
+ } else {
+ Sval = (*S1_ptr)(i, j, k);
+ }
+ int idxS = floor((Sval - S1_min) * hS_inv);
+ if (idxS < 0) {
+ idxS = 0;
+ } else if (idxS >= NS) {
+ idxS = 0;
+ }
+ local_hist[index(idxS, idxT, NS, NT)] += 1;
+ }
+ }
+ }
+
+ MPI_Barrier(MPI_COMM_WORLD); // Wait for everyone to finish
+ if (local_rank == 0) {
+ int *glob_hist = (int *)calloc(NS * NT, sizeof(int));
+ 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
+ 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;
+ }
+ }
+ free(glob_hist);
+ } else {
+ MPI_Reduce(local_hist, NULL, // send and receive buffers
+ NS * NT, MPI_INT, // count and datatype
+ MPI_SUM, 0, // Reduction operator and root process #
+ MPI_COMM_WORLD); // Communicator
+ }
+ free(local_hist);
+}
diff --git a/src/VERSION b/src/VERSION
index ae8b19475b67f0169e4f6057b5c5728ce674b540..81e4488aef22f8c969bdd8aada2b50890c93a412 100644
--- a/src/VERSION
+++ b/src/VERSION
@@ -1,3 +1,3 @@
MAJOR_VERSION=2
MINOR_VERSION=1
-PATCH_VERSION=8
+PATCH_VERSION=9
diff --git a/src/cases/derivatives/derivatives.cpp b/src/cases/derivatives/derivatives.cpp
index e34cf225c1f576bd8c7423d00db8723d79d3b4fb..99b158abab73d374d9a00df52328ea5a1cbe9dcd 100644
--- a/src/cases/derivatives/derivatives.cpp
+++ b/src/cases/derivatives/derivatives.cpp
@@ -8,6 +8,10 @@
#include "../../Options.hpp" // config-file parser
#include "../../Science.hpp" // Science content
+// Defines limits of intrinsic types. Used as default values for
+// T1_max, T1_min, S1_max and/orS1_min
+#include <limits>
+
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
@@ -29,6 +33,12 @@ 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;
+char T1_name, S1_name;
+int NT, NS;
+string QSP_filename;
+
// physical parameters
double visco; // viscosity (m^2/s)
double rho_0; // reference density (kg/m^3)
@@ -38,7 +48,7 @@ double g; // acceleration due to gravity (m/s^2)
int plotnum;
// Derivative options
-string deriv_filenames; // file name to take derivative of
+string deriv_filenames; // file name to take derivative of
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
@@ -52,8 +62,9 @@ bool do_enst_stretch; // Do enstrophy production term from vortex
bool do_Q; // Do second invariant of grad(u,v,w)?
bool do_R; // Do third invariant/det of grad(u,v,w)?
bool do_Q_and_R; // Do the second and third invariants?
-bool do_lambda2; // Do lambda2/Hussain's Lambda?
+bool do_lambda2; // Do lambda2/Hussain's Lambda?
bool v_exist; // Does the v field exist?
+bool do_hist; // Create QSP data?
/* ------------------ Adjust the class --------------------- */
@@ -87,7 +98,7 @@ class userControl : public BaseCase {
/* Set other things */
double get_visco() const { return visco; }
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) {
@@ -116,24 +127,25 @@ class userControl : public BaseCase {
//string prev_deriv, base_field;
vector<string> fields; // vector of fields to take derivatives
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 ) {
+ 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 )
+ if ( do_enst_stretch )
temp3 = alloc_array(Nx,Ny,Nz);
}
-
- if ( do_barvor ) {
+
+ if ( do_barvor ) {
temp4 = alloc_array(Nx,Ny,Nz);
}
if ( do_lambda2 ) {
- A11 = alloc_array(Nx,Ny,Nz);
- A12 = alloc_array(Nx,Ny,Nz);
- A13 = alloc_array(Nx,Ny,Nz);
- A22 = alloc_array(Nx,Ny,Nz);
- A23 = alloc_array(Nx,Ny,Nz);
- A33 = alloc_array(Nx,Ny,Nz);
+ A11 = alloc_array(Nx,Ny,Nz);
+ A12 = alloc_array(Nx,Ny,Nz);
+ A13 = alloc_array(Nx,Ny,Nz);
+ A22 = alloc_array(Nx,Ny,Nz);
+ A23 = alloc_array(Nx,Ny,Nz);
+ A33 = alloc_array(Nx,Ny,Nz);
}
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
@@ -154,7 +166,7 @@ class userControl : public BaseCase {
// parse for expansion type
find_expansion(grid_type, expan, fields[var_num]);
- // read the field and setup for derivative
+ // read the field and setup for derivative
if ( fields[var_num] == "v" ) {
saved_v = true;
init_tracer_restart(fields[var_num],v);
@@ -224,7 +236,7 @@ class userControl : public BaseCase {
}
}
}
-
+
// Baroclinic vorticity
if ( do_barvor ) {
@@ -239,23 +251,23 @@ class userControl : public BaseCase {
compute_baroclinic_vort_x(deriv_var, u, gradient_op, grid_type);
deriv_var=deriv_var*g;
write_array(deriv_var,"barx",plotnum);
- if (master()) fprintf(stdout,"Completed the write for barx.%d\n",plotnum);
+ if (master()) fprintf(stdout,"Completed the write for barx.%d\n",plotnum);
}
-
+
compute_baroclinic_vort_y(deriv_var, u, gradient_op, grid_type);
deriv_var=deriv_var*g;
write_array(deriv_var,"bary",plotnum);
if (master()) fprintf(stdout,"Completed the write for bary.%d\n",plotnum);
-
+
// Restart u
u = 0;
}
- // read in fields (if not already stored in memory)
+ // read in fields (if not already stored in memory)
if ( do_vor_x or do_vor_y or do_vor_z or
- do_enstrophy or do_dissipation or
- do_vort_stretch or do_enst_stretch or
- do_Q or do_R or do_Q_and_R or do_lambda2 ) {
+ do_enstrophy or do_dissipation or
+ do_vort_stretch or do_enst_stretch or
+ do_Q or do_R or do_Q_and_R or do_lambda2 or do_hist ) {
// u
init_tracer_restart("u",u);
// v
@@ -299,7 +311,7 @@ class userControl : public BaseCase {
write_array(deriv_var,"vortz",plotnum);
if (master())
fprintf(stdout,"Completed the write for vortz.%d\n",plotnum);
- }
+ }
// Calculate Enstrophy
if ( do_enstrophy ) {
@@ -373,7 +385,7 @@ class userControl : public BaseCase {
write_array(deriv_var,"Q",plotnum);
if (master()) fprintf(stdout,"Completed the write for Q.%d\n",plotnum);
}
-
+
// Calculate R/third invariant of grad(u,v,w)
if ( do_R or do_Q_and_R ) {
R_invt(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type, v_exist);
@@ -385,13 +397,29 @@ class userControl : public BaseCase {
// Calculate lambda2/second eigenvalue of S^2+Omega^2
if ( do_lambda2 ){
- compute_lambda2( deriv_var, u, v, w, *temp1, *temp2, gradient_op,
+ compute_lambda2( deriv_var, u, v, w, *temp1, *temp2, gradient_op,
grid_type, *A11, *A12, *A13, *A22, *A23, *A33);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max lam2: %.6g\n",max_var);
write_array(deriv_var,"lam2",plotnum);
if (master()) fprintf(stdout,"Completed the write for lam2.%d\n",plotnum);
-
+
+ }
+
+ // Compute QSP data. The code promises to not mutate the arrays,
+ // nor to make deep copies of them
+ if ( do_hist ){
+ // Read in T to temp1 if required.
+ if (T1_name == 't' || S1_name == 't' || T1_name == 'T' ||
+ 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 (master()) {
+ fprintf(stdout, "Completed the write for QSP.%d\n", plotnum);
+ }
}
}
}
@@ -465,6 +493,16 @@ 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("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_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_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);
@@ -489,13 +527,13 @@ int main(int argc, char ** argv) {
}
if (visco==1.0){
if (master())
- fprintf(stdout,"You may have forgotten to specify viscosity, "
+ fprintf(stdout,"You may have forgotten to specify viscosity, "
"Using default value visco = 1.\n");
}
if (rho_0==1000.0){
if (master())
- fprintf(stdout,"You may have forgotten to specify reference density, "
+ fprintf(stdout,"You may have forgotten to specify reference density, "
"Using default value rho_0 = 1000.\n");
}
diff --git a/src/cases/derivatives/spins.conf b/src/cases/derivatives/spins.conf
index 8b53d2d202b4df0f4f8e2952e8768967c15ba688..9e9304c8f827a3a96740c154769aa2182d8c29a5 100644
--- a/src/cases/derivatives/spins.conf
+++ b/src/cases/derivatives/spins.conf
@@ -42,3 +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
+QSP_filename = QSP
+NT = 10
+NS = 10