Commit c904026c authored by Christopher Subich's avatar Christopher Subich
Browse files

Merge branch 'master' into 'master'

Add QSPCount method to Science

See merge request SPINS/SPINS_main!12
parents d2801b45 d47b2ffa
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
......@@ -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
......
#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);
}
MAJOR_VERSION=2
MINOR_VERSION=1
PATCH_VERSION=8
PATCH_VERSION=9
......@@ -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,