Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • a2grace/SPINS_main
  • k2abduls/SPINS_main
  • klrowe/SPINS_main
  • dallum/SPINS_main
  • ncastrof/SPINS_main
  • spaceforce/SPINS_main
  • c2xu/SPINS_main
  • ddeepwel/SPINS_main
  • bastorer/SPINS_main
  • SPINS/SPINS_main
10 results
Show changes
Commits on Source (32)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1916 additions and 43 deletions
make_deps.sh
View file @ 48f54537
......@@ -70,12 +70,12 @@ if [ ! "$BUILD_BLITZ" = "yes" ]; then
else
echo "Building Blitz++"
# Download the Blitz tarball if necessary
if [ ! -e "blitz_2010.tgz" ]; then
mv redist_libs/blitz_2010.tgz ./
if [ ! -e "blitz_1.0.2.tar.gz" ]; then
mv redist_libs/blitz_1.0.2.tar.gz ./
fi
(tar -xzvf blitz_2010.tgz > /dev/null) || (echo "Untar of Blitz FAILED"; exit 1);
pushd blitz
(./configure --prefix="$CWD" --disable-fortran "${BLITZ_OPTIONS}" > /dev/null) && \
(tar -xzvf blitz_1.0.2.tar.gz > /dev/null) || (echo "Untar of Blitz FAILED"; exit 1);
pushd blitz-1.0.2
(autoreconf -vif && ./configure --prefix="$CWD" --disable-fortran "${BLITZ_OPTIONS}" > /dev/null) && \
(make lib > /dev/null) && \
pushd blitz && (make install > /dev/null) && popd && \
pushd lib && (make install > /dev/null) && popd && \
......@@ -92,14 +92,14 @@ if [ ! "$BUILD_FFTW" = "yes" ]; then
else
echo "Building FFTW"
# Download FFTW if necessary
if [ ! -e "fftw-3.3.3.tar.gz" ]; then
mv redist_libs/fftw-3.3.3.tar.gz ./
if [ ! -e "fftw-3.3.9.tar.gz" ]; then
mv redist_libs/fftw-3.3.9.tar.gz ./
fi
(tar -xzvf fftw-3.3.3.tar.gz > /dev/null)
(tar -xzvf fftw-3.3.9.tar.gz > /dev/null)
if [ 0 -ne $? ]; then
echo "Untar of FFTW FAILED"; exit 1
fi
pushd fftw-3.3.3
pushd fftw-3.3.9
# The "${FFTW_OPTIONS[@]}" syntax expands FFTW_OPTIONS as an array variable;
# this allows for multi-word arguments like 'CFLAGS="-O3 --fast-math"' to
# work properly as a single argument from configure's perspective.
......
matlab/spins_QSP_csv.m 0 → 100644
View file @ 48f54537
function [T1_max, T1_min, S1_max, S1_min, pdfCount] = spins_QSP_csv(filename)
%% Reads in the csv file and interprets it in the correct way for the user.
%%
%% Input:
%% Filename (string): name of csv file.
%%
%% Output:
%% T1_max (float): upper limit of the maximum-valued bin of T1 tracer
%% T1_min (float): lower limit of the minimum-valued bin of T1 tracer
%% S1_max (float): upper limit of the maximum-valued bin of S1 tracer
%% S1_min (float): lower limit of the minimum-valued bin of S1 tracer
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 = pdfCount / sum(sum(pdfCount));
end
src/ESolver.cpp
View file @ 48f54537
#include "TArray.hpp"
#include "blitz/array.h"
#include "blitz/tinyvec-et.h"
#include "ESolver.hpp"
#include "gmres_1d_solver.hpp"
#include "gmres_2d_solver.hpp"
......
src/Parformer.cpp
View file @ 48f54537
......@@ -2,7 +2,6 @@
#include "Parformer.hpp"
#include "Par_util.hpp"
#include <blitz/tinyvec-et.h>
#include <stdio.h>
#include <iostream>
......
src/Science.hpp
View file @ 48f54537
......@@ -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);
......@@ -130,9 +130,165 @@ void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type, bool v_exist);
inline double nleos_inline(double T, double S){
// 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) {
S = 0;
};
const double N0 = 9.99843699e+02;
const double N1 = 7.35212840e+00;
const double N2 = -5.45928211e-02;
const double N3 = 3.98476704e-04;
const double N4 = 2.96938239e+00;
const double N5 = -7.23268813e-03;
const double N6 = 2.12382341e-03;
const double N7 = 1.04004591e-02;
const double N8 = 1.03970529e-07;
const double N9 = 5.18761880e-06;
const double N10 = -3.24041825e-08;
const double N11 = -1.23869360e-11;
const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
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
const double D0 = 1.00000000e+00;
const double D1 = 7.28606739e-03;
const double D2 = -4.60835542e-05;
const double D3 = 3.68390573e-07;
const double D4 = 1.80809186e-10;
const double D5 = 2.14691708e-03;
const double D6 = -9.27062484e-06;
const double D7 = -1.78343643e-10;
const double D8 = 4.76534122e-06;
const double D9 = 1.63410736e-09;
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;
}
BZ_DECLARE_FUNCTION2(nleos_inline)
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
// 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;
}
BZ_DECLARE_FUNCTION(compute_alpha)
inline double compute_beta(double T0, double S0){
// 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;
double SpdS = S0 + dS;
double beta = (nleos_inline(T0,SpdS) - nleos_inline(T0,S0))/dS;
return beta;
}
BZ_DECLARE_FUNCTION(compute_beta)
inline double compute_rho0(double T0, double S0){
// Computes the reference density at S0 and T0
double rho_0 = nleos_inline(T0,S0);
return rho_0;
}
BZ_DECLARE_FUNCTION(compute_rho0)
void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
TArrayn::DTArray & S, const double & T0, const double & S0);
void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T);
void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0 = -10, double S0 = -2);
/*
inline double fresh_quad(double T){
// 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)
//0 pressure relative to sea surface pressure.
const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
const double C = -0.007641729398834; // Fitting Constant
return rho_max + C*(T - Tmd)*(T - Tmd);
}
BZ_DECLARE_FUNCTION(fresh_quad)
*/
// lambda2, second eigenvalue of S^2+Omega^2
void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
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);
/* 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).
*/
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)
//
inline double eqn_of_state(double T, double S){
// Returns the density anomaly (kg/m^3) for water at the given
// temperature T (degrees celsius) and salinity S (PSU?)
......@@ -163,7 +319,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
......
src/Science/QSPCount.cpp 0 → 100644
View file @ 48f54537
#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;
class QSPVector {
private:
double *data;
int length;
int rows;
int cols;
public:
explicit QSPVector(int length)
: data((double *)calloc(length, sizeof(double))), length(length), rows(0),
cols(0) {
if (!data) {
std::cout << "Error! Data could not be initialized.\n";
}
}
QSPVector(int Nx, int Ny)
: data((double *)calloc(Nx * Ny, sizeof(double))), length(Nx * Ny),
rows(Nx), cols(Ny) {
if (!data) {
std::cout << "Error! Data could not be initialized.\n";
}
}
double *raw() const { return data; }
int Length() const { return length; }
int Rows() const { return rows; }
int Cols() const { return cols; }
double operator[](int i) const {
assert(i >= 0 && i < length);
return data[i];
}
double operator()(int row, int col) const {
assert(row >= 0 && row < rows);
assert(col >= 0 && col < cols);
return data[(row * cols) + col];
}
double &operator()(int row, int col) {
assert(row >= 0 && row < rows);
assert(col >= 0 && col < cols);
return data[(row * cols) + col];
}
~QSPVector() { free(data); }
};
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
}
}
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;
}
TArrayn::DTArray *QSPPtr(const QSPData &qsp_data, const QSPType &type) {
TArrayn::DTArray *ptr = NULL;
switch (type) {
case QSP_u:
ptr = qsp_data.u;
break;
case QSP_v:
ptr = qsp_data.v;
break;
case QSP_w:
ptr = qsp_data.w;
break;
case QSP_ke: // This is an odd case.
break;
case QSP_rho:
ptr = qsp_data.rho;
break;
case QSP_temp:
ptr = qsp_data.temp;
break;
case QSP_salinity:
ptr = qsp_data.salinity;
break;
case QSP_custom: // Make sure to set the pointer yourself. Can't do it here.
break;
}
return ptr;
}
// 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();
// 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;
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));
}
}
void QSPCount(QSPOptions qsp_options, QSPData qsp_data) {
int local_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
// 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 =
(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(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 (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++) {
// min is set to the highest possible value so it always gets changed
double tmp_z_min = std::numeric_limits<double>::max();
// max is set to the lowest possible value so it always gets changed
double tmp_z_max = -tmp_z_min;
double tmp;
for (int kk = k_low; kk <= k_high; kk++) {
tmp = (*qsp_data.zgrid)(ii, jj, kk);
if (tmp > tmp_z_max) {
tmp_z_max = tmp;
} else if (tmp < tmp_z_min) {
tmp_z_min = tmp;
}
}
local_z_max(ii, jj) = tmp_z_max;
local_z_min(ii, jj) = tmp_z_min;
}
}
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
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++) {
switch (T1_type) {
case QSP_ke:
Tval = 0;
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;
break;
case QSP_rho:
tmp = (*T1_ptr)(i, j, k);
Tval = eqn_of_state_t(tmp);
break;
default:
Tval = (*T1_ptr)(i, j, k);
break;
}
switch (S1_type) {
case QSP_ke:
Sval = 0;
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;
break;
case QSP_rho:
tmp = (*S1_ptr)(i, j, k);
Sval = eqn_of_state_t(tmp);
break;
default:
Sval = (*S1_ptr)(i, j, k);
break;
}
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 (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 < 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 == 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)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
volume_weight = arc * (Lzmax_now - Lzmin_now);
} else {
volume_weight = 1.0;
}
local_hist(idxS, idxT) += volume_weight;
}
}
}
MPI_Barrier(MPI_COMM_WORLD); // Wait for everyone to finish
QSP_write(local_rank, local_hist, qsp_options, qsp_data);
}
src/Science/compute_lambda2.cpp 0 → 100644
View file @ 48f54537
#include "../Science.hpp"
#include <blitz/array.h>
#include <string>
#include <math.h>
using namespace Transformer;
using namespace blitz;
// Calculates the second largest eigenvalue of S^2+Omega^2
// S and Omega are the symmetric and antisymmetric components of u_i,j, respectively
// Note: S^2+Omega^2 = 0.5*(u_{i,k}u_{k,j}+u_{k,i}u_{j,k})
// This is used to find coherent structures correlated with eddies in turbulent flow
void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
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) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
TArrayn::DTArray * A_ref; //This will hold references to the components
TArrayn::DTArray * vel_mm; //This will hold the velocities
TArrayn::DTArray * vel_nn; //This will hold the velocities
std::string vel_labs[3] = {"u","v","w"};
int ctr = 0;
// Construct S^2+Omega^2
for (int mm=0; mm<3; mm++){
for (int nn=mm; nn<3; nn++){
// Choose which component will be constructed
if (ctr==0){A_ref = &A11;}
else if (ctr==1){A_ref = &A12;}
else if (ctr==2){A_ref = &A13;}
else if (ctr==3){A_ref = &A22;}
else if (ctr==4){A_ref = &A23;}
else if (ctr==5){A_ref = &A33;}
// A switch for the velocities
// Useful for computing the u_{i,k}u_{j,k} terms
if (mm==0){vel_mm = &u;}
else if (mm==1){vel_mm = &v;}
else if(mm==2){vel_mm = &w;}
if (nn==0){vel_nn = &u;}
else if (nn==1){vel_nn = &v;}
else if (nn==2){vel_nn = &w;}
// u_{i,k}u_{k,j}
// Get u_{i,x}
find_expansion(grid_type, expan, vel_labs[mm]);
gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp1,false);
// Get u_{x,j}
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
if (nn==0) { gradient_op->get_dx(&temp2,false); }
else if (nn==1) { gradient_op->get_dy(&temp2,false); }
else if (nn==2) { gradient_op->get_dz(&temp2,false); }
*A_ref = temp1 * temp2;
// Get u_{i,y}
find_expansion(grid_type, expan, vel_labs[mm]);
gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp1,false);
// Get u_{y,j}
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
if (nn==0) { gradient_op->get_dx(&temp2,false); }
else if (nn==1) { gradient_op->get_dy(&temp2,false); }
else if (nn==2) { gradient_op->get_dz(&temp2,false); }
*A_ref = temp1 * temp2;
// Get u_{i,z}
find_expansion(grid_type, expan, vel_labs[mm]);
gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get u_{z,j}
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
if (nn==0) { gradient_op->get_dx(&temp2,false); }
else if (nn==1) { gradient_op->get_dy(&temp2,false); }
else if (nn==2) { gradient_op->get_dz(&temp2,false); }
*A_ref = temp1 * temp2;
// u_{k,i}u_{j,k}
// Get u_{x,i}
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
if (mm==0) { gradient_op->get_dx(&temp1,false); }
else if (mm==1) { gradient_op->get_dy(&temp1,false); }
else if (mm==2) { gradient_op->get_dz(&temp1,false); }
// Get u_{j,x}
find_expansion(grid_type, expan, vel_labs[nn]);
gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
*A_ref = temp1 * temp2;
// Get u_{y,i}
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
if (mm==0) { gradient_op->get_dx(&temp1,false); }
else if (mm==1) { gradient_op->get_dy(&temp1,false); }
else if (mm==2) { gradient_op->get_dz(&temp1,false); }
// Get u_{j,y}
find_expansion(grid_type, expan, vel_labs[nn]);
gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
*A_ref = temp1 * temp2;
// Get u_{z,i}
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
if (mm==0) { gradient_op->get_dx(&temp1,false); }
else if (mm==1) { gradient_op->get_dy(&temp1,false); }
else if (mm==2) { gradient_op->get_dz(&temp1,false); }
// Get u_{j,z}
find_expansion(grid_type, expan, vel_labs[nn]);
gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
*A_ref = temp1 * temp2;
*A_ref = (*A_ref)*0.5;
ctr++;
}
}
// Now for the eigenvalue algorithm from:
// https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
// In the above article, we transform A into B=(1/p)*(A-qI)
// Then, lambda2 = q + 2*p*cos(arccos(det(B)/2)/3 + k*pi/3), k=0,1,2
// Define q
temp1 = (A11+A22+A33)/3.0;
// Define p
temp2 = pow2(A11-temp1) + pow2(A22-temp1) + pow2(A33-temp1) +2.0*(pow2(A12)+pow2(A13)+pow2(A23));
temp2 = sqrt(temp2/6.0);
// Transform A
A11 = (A11-temp1)/temp2;
A22 = (A22-temp1)/temp2;
A33 = (A33-temp1)/temp2;
A12 = A12/temp2;
A13 = A13/temp2;
A23 = A23/temp2;
// Calculate the determinant, using lambda2 as a dummy array
lambda2 = 0.5*(A11*(A22*A33-A23*A23)-A12*(A12*A33-A13*A23)+A13*(A12*A23-A13*A22));
// Since det(B)/2 feeds into acos, make sure it's between -1 and 1
for (int i = lambda2.lbound(blitz::firstDim); i <= lambda2.ubound(blitz::firstDim); i++) { // outer loop over slowest-varying dimension
for (int k = lambda2.lbound(blitz::thirdDim); k <= lambda2.ubound(blitz::thirdDim); k++) { // middle loop over next-slowest varying dimension
for (int j = lambda2.lbound(blitz::secondDim); j <= lambda2.ubound(blitz::secondDim); j++) { // inner loop over fastest varying dimeion
if (isnan(lambda2(i,j,k))) {lambda2(i,j,k) = 0;}
else if (lambda2(i,j,k) < -1.0) {lambda2(i,j,k) = -1.0;}
else if (lambda2(i,j,k) > 1.0) {lambda2(i,j,k) = 1.0;}
// first conditional works as a safety check because the det is NaN iff temp2=0; in this case the 2nd eig is temp1
}}}
// Now calculate the actual lamba2.
lambda2 = temp1+2.0*temp2*cos(acos(lambda2)/3.0+4.0*M_PI/3.0);
}
src/Science/eos.cpp 0 → 100644
View file @ 48f54537
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
#include <stdlib.h>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0, double S0) {
bool is_T0_valid = false;
bool is_S0_valid = false;
if (T0 >= -2) {
is_T0_valid = true;
}
if (S0 >= 0) {
is_S0_valid = true;
}
if (eos_type == "QUADEOS") {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the QUADEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call quadeos
quadeos(rho,T);
}
else if (eos_type == "LINEOS") {
if (is_T0_valid & is_S0_valid) {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the LINEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call lineos
lineos(rho,T,S,T0,S0);
}
else {
if (master()) {
fprintf(stderr,"Invalid option for background temperature or salinity. Exiting.\n");
}
MPI_Finalize(); exit(0);
}
}
else if (eos_type == "NLEOS") {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the NLEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call nleos
nleos(rho,T,S);
}
else {
if (master()) {
fprintf(stderr,"Invalid option received for eos type. Exiting.\n");
}
MPI_Finalize(); exit(0);
}
}
src/Science/lineos.cpp 0 → 100644
View file @ 48f54537
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S,
const double & T0, const double & S0) {
double alpha = compute_alpha(T0,S0);
double beta = compute_beta(T0,S0);
double rho_0 = compute_rho0(T0,S0);
rho = rho_0*(1 + alpha/rho_0*(T - T0) + beta/rho_0*(S - S0));
}
src/Science/nleos.cpp 0 → 100644
View file @ 48f54537
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S) {
// Returns the density (kg/m^3)
// MacDougall et. al. 2003 (JAOT 20)
//
// This EOS is a rational function taken at pressure = 0
//
// First find any regions of negative salinity and set them to zero for this operation.
//Numerator
const double N0 = 9.99843699e+02;
const double N1 = 7.35212840e+00;
const double N2 = -5.45928211e-02;
const double N3 = 3.98476704e-04;
const double N4 = 2.96938239e+00;
const double N5 = -7.23268813e-03;
const double N6 = 2.12382341e-03;
const double N7 = 1.04004591e-02;
const double N8 = 1.03970529e-07;
const double N9 = 5.18761880e-06;
const double N10 = -3.24041825e-08;
const double N11 = -1.23869360e-11;
const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
rho = 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
const double D0 = 1.00000000e+00;
const double D1 = 7.28606739e-03;
const double D2 = -4.60835542e-05;
const double D3 = 3.68390573e-07;
const double D4 = 1.80809186e-10;
const double D5 = 2.14691708e-03;
const double D6 = -9.27062484e-06;
const double D7 = -1.78343643e-10;
const double D8 = 4.76534122e-06;
const double D9 = 1.63410736e-09;
const double D10 = 5.30848875e-06;
const double D11 = -3.03175128e-16;
const double D12 = -1.27934137e-17;
rho = rho/(D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(max(S,0))*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p)));
}
src/Science/quadeos.cpp 0 → 100644
View file @ 48f54537
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T) {
// 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)
//0 pressure relative to sea surface pressure.
const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
const double C = -0.007641729398834; // Fitting Constant
rho = rho_max + C*(T - Tmd)*(T - Tmd);
}
src/Sorter.cpp
View file @ 48f54537
......@@ -10,7 +10,6 @@
#include <vector>
#include <algorithm>
#include <stdlib.h>
#include <blitz/tinyvec-et.h>
#include "Sorter.hpp"
......
src/Splits.cpp
View file @ 48f54537
......@@ -2,7 +2,6 @@
#include <vector>
#include <blitz/array.h>
#include <blitz/tinyvec-et.h>
#include <stdio.h>
#include <iostream>
#include <complex>
......@@ -121,7 +120,7 @@ src_temp(0), dst_temp(0), issue_warning(true) {
rec_types.resize(num_proc); // make sure our vector can hold everything
for (int i = 0; i < num_proc; i++) {
// We'll have to use the MPI_Type_struct syntax, so build arrays:
/*// We'll have to use the MPI_Type_struct syntax, so build arrays:
MPI_Datatype types[2] = {vec_type,MPI_UB};
int counts[2]; counts[0] = extent_from[i]; counts[1] = 1;
// Now, set displacements to set the "end" of the datatype a full
......@@ -129,8 +128,23 @@ src_temp(0), dst_temp(0), issue_warning(true) {
MPI_Aint displs[2] = {0,sizeof(T)*sizes[untouched_dim]*sizes[from_dim]};
// And make the type
MPI_Type_struct(2,counts,displs,types,&rec_types[i]);
MPI_Type_struct(2,counts,displs,types,&rec_types[i]);*/
// This code formerly used the MPI_UB marker to define the upper bound
// of the derived datatype. This was deprecated as part of the MPI-2
// standard, and it was removed in MPI-3. The replacement is to use
// MPI_Type_create_resized to create a resized type directly.
// First, create the multi-vector containing the data:
MPI_Datatype multivec;
MPI_Type_contiguous(extent_from[i],vec_type,&multivec);
MPI_Type_create_resized(multivec, 0, // base type, lower bound
sizeof(T)*sizes[untouched_dim]*sizes[from_dim], // extent
&rec_types[i]); // new type
MPI_Type_commit(&rec_types[i]);
MPI_Type_free(&multivec);
/* Impelementation note: one optimization possible here is to collapse
this typelist. This impelementation makes no assumptions about how
many part-rows are received per processor, but our default splitting
......
src/TArray.cpp
View file @ 48f54537
#include "TArray.hpp"
#include "Plan_holder.hpp"
#include <blitz/tinyvec-et.h>
#include <assert.h>
#include <fftw3.h> // FFTW
/* Implementation of TArray transform functions for the relevant cases.
......
src/T_util.cpp
View file @ 48f54537
#include "T_util.hpp"
#include "TArray.hpp"
#include <blitz/array.h>
#include <blitz/tinyvec-et.h>
#include "Par_util.hpp"
#include <complex>
#include <string>
......
src/VERSION
View file @ 48f54537
MAJOR_VERSION=2
MINOR_VERSION=1
PATCH_VERSION=6
PATCH_VERSION=9
src/cases/brydon_test/brydon_test.cpp 0 → 100644
View file @ 48f54537
This diff is collapsed.
src/cases/brydon_test/spins.conf 0 → 100644
View file @ 48f54537
## Salt-temperature configuration file
# Spatial Parameters
Lx = 1.0
Ly = 1.0
Lz = 0.1
Nx = 1024
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_T = 2e-6
kappa_S = 1e-7
# Problem Parameters (delta_rho is percentage)
delta_T = 25.0
delta_S = 5.0
T_0 = 15.0
S_0 = 0.0
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 100
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 1e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
compute_BPE = false
compute_internal_to_BPE = false
#compute_stresses_top = false
#compute_stresses_bottom = false
src/cases/derivatives/derivatives.cpp
View file @ 48f54537
......@@ -8,6 +8,11 @@
#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 <cassert>
#include <limits>
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
......@@ -29,6 +34,13 @@ 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;
string T1_name, S1_name;
int NT, NS;
string QSP_filename;
bool use_salinity;
// physical parameters
double visco; // viscosity (m^2/s)
double rho_0; // reference density (kg/m^3)
......@@ -38,7 +50,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,7 +64,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 v_exist; // Does the v field exist?
bool do_hist; // Create QSP data?
/* ------------------ Adjust the class --------------------- */
......@@ -63,6 +77,8 @@ class userControl : public BaseCase {
DTArray deriv_var; // array for derivative
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; }
......@@ -85,7 +101,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) {
......@@ -114,17 +130,43 @@ 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 ) {
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 (mapped) {
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 or do_hist ) 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);
}
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
plotnum = plotnum + step_sequence ) {
......@@ -144,7 +186,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);
......@@ -214,7 +256,7 @@ class userControl : public BaseCase {
}
}
}
// Baroclinic vorticity
if ( do_barvor ) {
......@@ -229,23 +271,23 @@ class userControl : public BaseCase {
compute_baroclinic_vort_x(deriv_var, u, gradient_op, grid_type);
deriv_var=deriv_var*g/rho_0;
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/rho_0;
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 ) {
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
......@@ -289,7 +331,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 ) {
......@@ -363,7 +405,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);
......@@ -372,6 +414,82 @@ class userControl : public BaseCase {
write_array(deriv_var,"R",plotnum);
if (master()) fprintf(stdout,"Completed the write for R.%d\n",plotnum);
}
// Calculate lambda2/second eigenvalue of S^2+Omega^2
if ( do_lambda2 ){
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 ){
// If user asked to grab the grids, we populate the grids
// with the correct data from disk
if (mapped) {
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);
}
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);
}
}
}
}
......@@ -419,7 +537,7 @@ int main(int argc, char ** argv) {
add_option("type_z",&zgrid_type,"Grid type in Z");
option_category("Physical parameters");
add_option("visco",&visco,0.0,"Viscosity");
add_option("visco",&visco,1.0,"Viscosity");
add_option("rho_0",&rho_0,1000.0,"Reference Density");
add_option("g",&g,9.81,"Acceleration due to gravity");
......@@ -443,6 +561,18 @@ int main(int argc, char ** argv) {
add_option("do_Q",&do_Q,false,"Calculate Q?");
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,"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("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");
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);
......@@ -465,6 +595,17 @@ int main(int argc, char ** argv) {
fprintf(stdout,"Simulation is 2 dimensional, "
"Ly has been changed to 1.0 for normalization.\n");
}
if (visco==1.0){
if (master())
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, "
"Using default value rho_0 = 1000.\n");
}
/* ------------------ Do stuff --------------------- */
userControl mycode; // Create an instantiated object of the above class
......
src/cases/derivatives/spins.conf
View file @ 48f54537
......@@ -40,4 +40,13 @@ do_enst_stretch = false
do_barvor = false
do_Q = false
do_R = false
do_Q_and_R = true
do_Q_and_R = false
do_lambda2 = false
# QSP Parameters
do_hist = false
T1 = temp
S1 = ke
QSP_filename = QSP
NT = 10
NS = 10