diff --git a/make_deps.sh b/make_deps.sh
index 5c8da679b73ad9407fc7e76aa1bd96a713aac737..c1c03a82873de41bb88ddbe4081adbac53b24e59 100755
--- a/make_deps.sh
+++ b/make_deps.sh
@@ -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.
diff --git a/matlab/spins_QSP_csv.m b/matlab/spins_QSP_csv.m
new file mode 100644
index 0000000000000000000000000000000000000000..d544e59fa3bae91a8c6c2480384a0ae9f6cfae0b
--- /dev/null
+++ b/matlab/spins_QSP_csv.m
@@ -0,0 +1,19 @@
+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
diff --git a/src/ESolver.cpp b/src/ESolver.cpp
index 656e4ea7e9db368bedddd3f518036739a538939f..4eab40950d4ee123d7365c986312ec4469d1aa1d 100644
--- a/src/ESolver.cpp
+++ b/src/ESolver.cpp
@@ -1,6 +1,5 @@
#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"
diff --git a/src/Parformer.cpp b/src/Parformer.cpp
index 331a2922a7814e9f108461a76bb33d7c5558f6b5..383bde710afa35632db5c47b0a31495e5543eba2 100644
--- a/src/Parformer.cpp
+++ b/src/Parformer.cpp
@@ -2,7 +2,6 @@
#include "Parformer.hpp"
#include "Par_util.hpp"
-#include <blitz/tinyvec-et.h>
#include <stdio.h>
#include <iostream>
diff --git a/src/Science.hpp b/src/Science.hpp
index 569b9a6c6c494ddbb583efef4e924702b6e43237..cd5c38f754ead297aca9c7e09f6a338d2d164096 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);
@@ -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
diff --git a/src/Science/QSPCount.cpp b/src/Science/QSPCount.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25fdae08abca743c8d0c48e222fd7722e134f65a
--- /dev/null
+++ b/src/Science/QSPCount.cpp
@@ -0,0 +1,458 @@
+#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);
+}
diff --git a/src/Science/compute_lambda2.cpp b/src/Science/compute_lambda2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abd08d62d940c41c535bfe9a72b3e48cc179ff42
--- /dev/null
+++ b/src/Science/compute_lambda2.cpp
@@ -0,0 +1,161 @@
+#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);
+
+
+}
diff --git a/src/Science/eos.cpp b/src/Science/eos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..77010f023921bc0cc23f1a8e76a4bace88aae211
--- /dev/null
+++ b/src/Science/eos.cpp
@@ -0,0 +1,84 @@
+#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);
+ }
+}
+
+
diff --git a/src/Science/lineos.cpp b/src/Science/lineos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fe65663a4eef4f0a39a35e98d6118e3327854fb1
--- /dev/null
+++ b/src/Science/lineos.cpp
@@ -0,0 +1,26 @@
+#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));
+
+}
diff --git a/src/Science/nleos.cpp b/src/Science/nleos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e270b9e571fbe956767c7afaa6b3e78af9a644ce
--- /dev/null
+++ b/src/Science/nleos.cpp
@@ -0,0 +1,63 @@
+#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)));
+}
+
diff --git a/src/Science/quadeos.cpp b/src/Science/quadeos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e5fb8e17b295fbaa8dc4695c7b353edd1cdce647
--- /dev/null
+++ b/src/Science/quadeos.cpp
@@ -0,0 +1,30 @@
+#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);
+
+}
diff --git a/src/Sorter.cpp b/src/Sorter.cpp
index 27836927fd8cf6304b09a56fe88c988b0ad1108e..21b008b84e61bc5e0bab79f1f086418f81056639 100644
--- a/src/Sorter.cpp
+++ b/src/Sorter.cpp
@@ -10,7 +10,6 @@
#include <vector>
#include <algorithm>
#include <stdlib.h>
-#include <blitz/tinyvec-et.h>
#include "Sorter.hpp"
diff --git a/src/Splits.cpp b/src/Splits.cpp
index 1e366bee7578ab2e7cd6f5f8fa87ff0abf8240b7..afb5162e76e7404e68b999170432d4e713c305e2 100644
--- a/src/Splits.cpp
+++ b/src/Splits.cpp
@@ -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
diff --git a/src/TArray.cpp b/src/TArray.cpp
index 5137e9246470d394ba6f4945bc1e9d2906848bb3..d10b09d13d04ee926293b0810251afb4115d04ff 100644
--- a/src/TArray.cpp
+++ b/src/TArray.cpp
@@ -1,6 +1,5 @@
#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.
diff --git a/src/T_util.cpp b/src/T_util.cpp
index 7503f132e5c428f1f71124a334c87f4fa61a8e67..7913ee76aafaddfba21dc5224e217ab9f81b0ab0 100644
--- a/src/T_util.cpp
+++ b/src/T_util.cpp
@@ -1,7 +1,6 @@
#include "T_util.hpp"
#include "TArray.hpp"
#include <blitz/array.h>
-#include <blitz/tinyvec-et.h>
#include "Par_util.hpp"
#include <complex>
#include <string>
diff --git a/src/VERSION b/src/VERSION
index 2944c9e8063acc3ac0de00f2350a9e70fc6318e9..81e4488aef22f8c969bdd8aada2b50890c93a412 100644
--- a/src/VERSION
+++ b/src/VERSION
@@ -1,3 +1,3 @@
MAJOR_VERSION=2
MINOR_VERSION=1
-PATCH_VERSION=6
+PATCH_VERSION=9
diff --git a/src/cases/brydon_test/brydon_test.cpp b/src/cases/brydon_test/brydon_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fae64e442dca599650b504ddd155402e800e91ad
--- /dev/null
+++ b/src/cases/brydon_test/brydon_test.cpp
@@ -0,0 +1,650 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+ w_f = -g*(eqn_of_state(*tracers[TEMP],*tracers[SALT]))/rho_0;
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/brydon_test/spins.conf b/src/cases/brydon_test/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c278245ccd17011b34d7b7cd6cd4dbdcef1ff
--- /dev/null
+++ b/src/cases/brydon_test/spins.conf
@@ -0,0 +1,67 @@
+## 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
diff --git a/src/cases/derivatives/derivatives.cpp b/src/cases/derivatives/derivatives.cpp
index ab383f08ef2f6d100e50a985837abb1b193969a9..3ed017656a9a31b9b916fc460af83303bc6ec5eb 100644
--- a/src/cases/derivatives/derivatives.cpp
+++ b/src/cases/derivatives/derivatives.cpp
@@ -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
diff --git a/src/cases/derivatives/spins.conf b/src/cases/derivatives/spins.conf
index 352e07a7d058b22f19f66eb1deb0da4c008cb935..3d05b5ba8bed408714ff6da43f738ac245dee918 100644
--- a/src/cases/derivatives/spins.conf
+++ b/src/cases/derivatives/spins.conf
@@ -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
diff --git a/src/cases/jmd_test/jmd_test.cpp b/src/cases/jmd_test/jmd_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..99ecd95bf751660bb07d77df3c9ab730f76e8a7b
--- /dev/null
+++ b/src/cases/jmd_test/jmd_test.cpp
@@ -0,0 +1,650 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+ w_f = -g*(nleos(*tracers[TEMP],*tracers[SALT]) - rho_0)/rho_0;
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/jmd_test/spins.conf b/src/cases/jmd_test/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c278245ccd17011b34d7b7cd6cd4dbdcef1ff
--- /dev/null
+++ b/src/cases/jmd_test/spins.conf
@@ -0,0 +1,67 @@
+## 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
diff --git a/src/cases/lineos_test/lineos_test.cpp b/src/cases/lineos_test/lineos_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80d7075250928ade6e372f1362341d58995da3ea
--- /dev/null
+++ b/src/cases/lineos_test/lineos_test.cpp
@@ -0,0 +1,655 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+
+ //Do Eos stuff
+ lineos(*temp1,*tracers[TEMP],*tracers[SALT],T_0,S_0);
+ rho_0 = compute_rho0(T_0,S_0);
+
+ w_f = -g/rho_0*(*temp1-rho_0);
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/lineos_test/spins.conf b/src/cases/lineos_test/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c278245ccd17011b34d7b7cd6cd4dbdcef1ff
--- /dev/null
+++ b/src/cases/lineos_test/spins.conf
@@ -0,0 +1,67 @@
+## 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
diff --git a/src/cases/nleos_test/nleos_test.cpp b/src/cases/nleos_test/nleos_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6bea267ca530f3bd3d07ede918adeeb533714faf
--- /dev/null
+++ b/src/cases/nleos_test/nleos_test.cpp
@@ -0,0 +1,655 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+
+ //Do Eos stuff
+ nleos(*temp1,*tracers[TEMP],*tracers[SALT]);
+ rho_0 = compute_rho0(T_0,S_0);
+
+ w_f = -g/rho_0*(*temp1-rho_0);
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/nleos_test/spins.conf b/src/cases/nleos_test/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c278245ccd17011b34d7b7cd6cd4dbdcef1ff
--- /dev/null
+++ b/src/cases/nleos_test/spins.conf
@@ -0,0 +1,67 @@
+## 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
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
diff --git a/src/cases/quadeos_test/nleos_test.cpp b/src/cases/quadeos_test/nleos_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6bea267ca530f3bd3d07ede918adeeb533714faf
--- /dev/null
+++ b/src/cases/quadeos_test/nleos_test.cpp
@@ -0,0 +1,655 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+
+ //Do Eos stuff
+ nleos(*temp1,*tracers[TEMP],*tracers[SALT]);
+ rho_0 = compute_rho0(T_0,S_0);
+
+ w_f = -g/rho_0*(*temp1-rho_0);
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/quadeos_test/quadeos_test.cpp b/src/cases/quadeos_test/quadeos_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d70c550917b18f3268b1f123eea62110d01f636a
--- /dev/null
+++ b/src/cases/quadeos_test/quadeos_test.cpp
@@ -0,0 +1,656 @@
+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
+#include "../../Options.hpp" // config-file parser
+#include <random/normal.h> // Blitz random number generator
+
+using namespace ranlib;
+
+// 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
+double MinX, MinY, MinZ; // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco; // viscosity (m^2/s)
+double mu; // dynamic viscosity (kg/(m·s))
+double kappa_T; // diffusivity of temperature (m^2/s)
+double kappa_S; // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2; // number of tracers (temperature, salt and dyes)
+const int TEMP = 0; // index for temperature
+const int SALT = 1; // index for salt
+
+// Problem parameters
+double delta_T; // temperature difference between different layers
+double delta_S; // salinity difference between different layers
+double T_0; // background temperature
+double S_0; // background salinity
+double delta_x; // horizontal transition length (m)
+double Lmix; // Width of mixed region (m)
+
+// Topography parameters
+double hill_height; // height of hill (m)
+double hill_centre; // position of hill peak (m)
+double hill_width; // width of hill (m)
+
+// Temporal parameters
+double final_time; // Final time (s)
+double plot_interval; // Time between field writes (s)
+double dt_max; // maximum time step (s)
+
+// Restarting options
+bool restarting; // are you restarting?
+double initial_time; // initial start time of simulation
+int restart_sequence; // output number to restart from
+
+// Dump parameters
+bool restart_from_dump; // restarting from dump?
+double compute_time; // requested computation time
+double avg_write_time; // average time to write all output fields at one output
+double real_start_time; // real (clock) time when simulation begins
+double compute_start_time; // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb; // Initial velocity perturbation
+bool compute_enstrophy; // Compute enstrophy?
+bool compute_dissipation; // Compute dissipation?
+bool compute_BPE; // Compute background potential energy?
+bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
+bool compute_stresses_top; // Compute top surface stresses?
+bool compute_stresses_bottom; // Compute bottom surface stresses?
+bool write_pressure; // Write the pressure field?
+int iter = 0; // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+ public:
+ // Grid and topography arrays
+ DTArray *zgrid; // Full grid fields
+ Array<double,1> xx, yy, zz; // 1D grid vectors
+ Array<double,1> topo; // topography vector
+ DTArray *Hprime; // derivative of topography vector
+
+ // Arrays and operators for derivatives
+ Grad * gradient_op;
+ DTArray *temp1, *dxdydz;
+
+ // Timing variables (for outputs and measuring time steps)
+ int plot_number; // plot output number
+ double next_plot; // time of next output write
+ double comp_duration; // clock time since computation began
+ double clock_time; // current clock time
+
+ // 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; }
+
+ // Record the gradient-taking object
+ void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+ // Coriolis parameter, viscosity, and diffusivities
+ double get_rot_f() const { return rot_f; }
+ double get_visco() const { return visco; }
+ double get_diffusivity(int t_num) const {
+ switch (t_num) {
+ case TEMP:
+ return kappa_T;
+ case SALT:
+ return kappa_S;
+ default:
+ if (master()) fprintf(stderr,"Invalid tracer number!\n");
+ MPI_Finalize(); exit(1);
+ }
+ }
+
+ // Temporal parameters
+ double init_time() const { return initial_time; }
+ int get_restart_sequence() const { return restart_sequence; }
+ double get_dt_max() const { return dt_max; }
+ double get_next_plot() { return next_plot; }
+
+ // Number of tracers (the first is density)
+ int numtracers() const { return Num_tracers; }
+
+ // Create mapped grid
+ bool is_mapped() const { return mapped; }
+ void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+ zgrid = alloc_array(Nx,Ny,Nz);
+
+ // over-write zz to be between -1 and 1
+ // (zz defined in automatic_grid below)
+ zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical
+
+ // Define topography
+ topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+ // make full grids
+ xg = xx(ii) + 0*jj + 0*kk;
+ yg = yy(jj) + 0*ii + 0*kk;
+ zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+ *zgrid = zg;
+
+ // Write the arrays and matlab readers
+ write_array(xg,"xgrid");
+ write_reader(xg,"xgrid",false);
+ if (Ny > 1) {
+ write_array(yg,"ygrid");
+ write_reader(yg,"ygrid",false);
+ }
+ write_array(zg,"zgrid");
+ write_reader(zg,"zgrid",false);
+ }
+
+ /* Initialize velocities */
+ void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+ if (master()) fprintf(stdout,"Initializing velocities\n");
+ // if restarting
+ if (restarting and !restart_from_dump) {
+ init_vels_restart(u, v, w);
+ } else if (restarting and restart_from_dump) {
+ init_vels_dump(u, v, w);
+ } else{
+ // else have a near motionless field
+ u = 0;
+ v = 0;
+ w = 0;
+ // Add a random perturbation to trigger any 3D instabilities
+ int myrank;
+ MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+ Normal<double> rnd(0,1);
+ for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+ rnd.seed(i);
+ for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+ for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+ u(i,j,k) += perturb*rnd.random();
+ w(i,j,k) += perturb*rnd.random();
+ if (Ny > 1 || rot_f != 0)
+ v(i,j,k) += perturb*rnd.random();
+ }
+ }
+ }
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0) {
+ write_array(v,"v",plot_number);
+ }
+ }
+ }
+
+ /* Initialize the tracers (density and dyes) */
+ void init_tracers(vector<DTArray *> & tracers) {
+ if (master()) fprintf(stdout,"Initializing tracers\n");
+ DTArray & temp = *tracers[TEMP];
+ DTArray & salt = *tracers[SALT];
+
+ if (restarting and !restart_from_dump) {
+ init_tracer_restart("t",temp);
+ init_tracer_restart("s",salt);
+ } else if (restarting and restart_from_dump) {
+ init_tracer_dump("t",temp);
+ init_tracer_dump("s",salt);
+ } else {
+ // temperature configuration
+ temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // salt configuration
+ salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+ // Important! if mapped, and t or s depends on z
+ // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+ // Write the array
+ write_array(temp,"t",plot_number);
+ write_array(salt,"s",plot_number);
+ }
+ }
+
+ /* Forcing in the momentum equations */
+ void forcing(double t, const DTArray & u, DTArray & u_f,
+ const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+ vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+ u_f = +rot_f*v;
+ v_f = -rot_f*u;
+
+ //Do Eos stuff
+ quadeos(*temp1,*tracers[TEMP]);
+ //The line below is patchwork, but is enough for a test case
+ rho_0 = compute_rho0(3.98,0);
+
+ w_f = -g/rho_0*(*temp1-rho_0);
+ *tracers_f[TEMP] = 0;
+ *tracers_f[SALT] = 0;
+ }
+
+ /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+ void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers, DTArray & pressure) {
+ // Set-up
+ if ( iter == 0 ) {
+ if ( compute_enstrophy or compute_dissipation or
+ compute_stresses_top or compute_stresses_bottom ) {
+ temp1 = alloc_array(Nx,Ny,Nz);
+ }
+ if ( compute_stresses_top or compute_stresses_bottom ) {
+ // initialize the vector of the bottom slope (Hprime)
+ Hprime = alloc_array(Nx,Ny,1);
+ if (is_mapped()) {
+ bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+ } else {
+ topo = 0*ii;
+ *Hprime = 0*ii + 0*jj;
+ }
+ }
+ // Determine last plot if restarting from the dump file
+ if (restart_from_dump) {
+ next_plot = (restart_sequence+1)*plot_interval;
+ }
+ // initialize the size of each voxel
+ dxdydz = alloc_array(Nx,Ny,Nz);
+ *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+ if (is_mapped()) {
+ *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+ }
+ }
+ // update clocks
+ if (master()) {
+ clock_time = MPI_Wtime();
+ comp_duration = clock_time - compute_start_time;
+ }
+
+ /* Calculate and write out useful information */
+
+ // Energy (PE assumes density is density anomaly)
+ double ke_x = 0, ke_y = 0, ke_z = 0;
+ if ( Nx > 1 ) {
+ ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+ }
+ if (Ny > 1 || rot_f != 0) {
+ ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+ }
+ if ( Nz > 1 ) {
+ ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+ }
+ double pe_tot;
+ if (is_mapped()) {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+ } else {
+ pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+ }
+ double BPE_tot = 0;
+ if (compute_BPE) {
+ // full energy budget for salt/heat stratifcations is incomplete
+ //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+ // g, rho_0, iter, true, is_mapped(), topo);
+ }
+ // Conversion from internal energy to background potential energy
+ double phi_i = 0;
+ if (compute_internal_to_BPE) {
+ // this is not finished for temperature/salt stratified cases
+ //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+ }
+ // viscous dissipation
+ double diss_tot = 0;
+ double max_diss = 0;
+ if (compute_dissipation) {
+ dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+ max_diss = psmax(max(*temp1));
+ diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+ }
+ // Vorticity / Enstrophy
+ double max_vort_x = 0, enst_x_tot = 0;
+ double max_vort_y = 0, enst_y_tot = 0;
+ double max_vort_z = 0, enst_z_tot = 0;
+ if (compute_enstrophy) {
+ // x-vorticity
+ if (Ny > 1 and Nz > 1) {
+ compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+ max_vort_x = psmax(max(abs(*temp1)));
+ enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // y-vorticity
+ if (Nx > 1 and Nz > 1) {
+ compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+ max_vort_y = psmax(max(abs(*temp1)));
+ enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ // z-vorticity
+ if (Nx > 1 and Ny > 1) {
+ compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+ max_vort_z = psmax(max(abs(*temp1)));
+ enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+ }
+ }
+ // max of fields
+ double max_u = psmax(max(abs(u)));
+ double max_v = psmax(max(abs(v)));
+ double max_w = psmax(max(abs(w)));
+ double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+ double max_temp = psmax(max(*tracers[TEMP]));
+ double min_temp = psmin(min(*tracers[TEMP]));
+ double max_salt = psmax(max(*tracers[SALT]));
+ double min_salt = psmin(min(*tracers[SALT]));
+ double max_rho = psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ double min_rho = psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+ // total mass (tracers[RHO] is non-dimensional density)
+ double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+ if (master()) {
+ // add diagnostics to buffers
+ string header, line;
+ add_diagnostic("Iter", iter, header, line);
+ add_diagnostic("Clock_time", comp_duration, header, line);
+ add_diagnostic("Time", time, header, line);
+ add_diagnostic("Max_vel", max_vel, header, line);
+ add_diagnostic("Max_temperature", max_temp, header, line);
+ add_diagnostic("Min_temperature", min_temp, header, line);
+ add_diagnostic("Max_salinity", max_salt, header, line);
+ add_diagnostic("Min_salinity", min_salt, header, line);
+ add_diagnostic("Max_density", max_rho, header, line);
+ add_diagnostic("Min_density", min_rho, header, line);
+ add_diagnostic("Mass", mass, header, line);
+ add_diagnostic("PE_tot", pe_tot, header, line);
+ if (compute_BPE) {
+ add_diagnostic("BPE_tot", BPE_tot, header, line);
+ }
+ if (compute_internal_to_BPE) {
+ add_diagnostic("BPE_from_int", phi_i, header, line);
+ }
+ if (compute_dissipation) {
+ add_diagnostic("Max_diss", max_diss, header, line);
+ add_diagnostic("Diss_tot", diss_tot, header, line);
+ }
+ if (Nx > 1) {
+ add_diagnostic("Max_u", max_u, header, line);
+ add_diagnostic("KE_x", ke_x, header, line);
+ }
+ if (Ny > 1 || rot_f != 0) {
+ add_diagnostic("Max_v", max_v, header, line);
+ add_diagnostic("KE_y", ke_y, header, line);
+ }
+ if (Nz > 1) {
+ add_diagnostic("Max_w", max_w, header, line);
+ add_diagnostic("KE_z", ke_z, header, line);
+ }
+ if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+ add_diagnostic("Max_vort_x", max_vort_x, header, line);
+ }
+ if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+ add_diagnostic("Max_vort_y", max_vort_y, header, line);
+ }
+ if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+ add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+ add_diagnostic("Max_vort_z", max_vort_z, header, line);
+ }
+
+ // Write to file
+ if (!(restarting and iter==0))
+ write_diagnostics(header, line, iter, restarting);
+ // and to the log file
+ fprintf(stdout,"[%d] (%.4g) %.4f: "
+ "%.4g %.4g %.4g %.4g\n",
+ iter,comp_duration,time,
+ max_u,max_v,max_w,max_rho);
+ }
+
+ // Top Surface Stresses
+ if ( compute_stresses_top ) {
+ stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+ // Bottom Surface Stresses
+ if ( compute_stresses_bottom ) {
+ stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+ }
+
+ /* Write to disk if at correct time */
+ if ((time - next_plot) > -1e-6) {
+ plot_number++;
+ comp_duration = MPI_Wtime(); // time just before write (for dump)
+ // Write the arrays
+ write_array(u,"u",plot_number);
+ write_array(w,"w",plot_number);
+ if (Ny > 1 || rot_f != 0)
+ write_array(v,"v",plot_number);
+ write_array(*tracers[TEMP],"t",plot_number);
+ write_array(*tracers[SALT],"s",plot_number);
+ if (write_pressure)
+ write_array(pressure,"p",plot_number);
+ // update next plot time
+ next_plot = next_plot + plot_interval;
+
+ // Find average time to write (for dump)
+ clock_time = MPI_Wtime(); // time just after write
+ avg_write_time = (avg_write_time*(plot_number-restart_sequence-1)
+ + (clock_time - comp_duration))/(plot_number-restart_sequence);
+ // Print information about plot outputs
+ write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+ }
+
+ // see if close to end of compute time and dump
+ check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+ plot_number, iter, u, v, w, tracers);
+ // Change dump log file if successfully reached final time
+ successful_dump(plot_number, final_time, plot_interval);
+ // increase counter
+ iter++;
+ }
+
+ // User specified variables to dump
+ void write_variables(DTArray & u,DTArray & v, DTArray & w,
+ vector<DTArray *> & tracers) {
+ write_array(u,"u.dump");
+ write_array(v,"v.dump");
+ write_array(w,"w.dump");
+ write_array(*tracers[TEMP],"t.dump");
+ write_array(*tracers[SALT],"s.dump");
+ }
+
+ // Constructor: Initialize local variables
+ userControl():
+ xx(split_range(Nx)), yy(Ny), zz(Nz),
+ topo(split_range(Nx)), gradient_op(0),
+ plot_number(restart_sequence),
+ next_plot(initial_time + plot_interval)
+ { compute_quadweights(
+ size_x(), size_y(), size_z(),
+ length_x(), length_y(), length_z(),
+ type_x(), type_y(), type_z());
+ // Create one-dimensional arrays for the coordinates
+ automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+ }
+};
+
+/* 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);
+
+ real_start_time = MPI_Wtime(); // start of simulation (for dump)
+ /* ------------------ 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");
+ add_option("min_x",&MinX,0.0,"Minimum X-value");
+ add_option("min_y",&MinY,0.0,"Minimum Y-value");
+ add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+ option_category("Grid expansion options");
+ 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");
+ add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+ option_category("Topography parameters");
+ add_option("hill_height",&hill_height,"Height of hill");
+ add_option("hill_centre",&hill_centre,"location of hill peak");
+ add_option("hill_width",&hill_width,"Width of hill");
+
+ option_category("Physical parameters");
+ add_option("g",&g,9.81,"Gravitational acceleration");
+ add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+ add_option("rho_0",&rho_0,1000.0,"Reference density");
+ add_option("visco",&visco,"Viscosity");
+ add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+ add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+ option_category("Problem parameters");
+ add_option("delta_T",&delta_T,"Temperature difference");
+ add_option("delta_S",&delta_S,"Salinity difference");
+ add_option("T_0",&T_0,"Background temperature");
+ add_option("S_0",&S_0,"Background salinity");
+ add_option("delta_x",&delta_x,"Horizontal transition half-width");
+ add_option("Lmix",&Lmix,"Width of collapse region");
+
+ option_category("Temporal options");
+ add_option("final_time",&final_time,"Final time");
+ add_option("plot_interval",&plot_interval,"Time between writes");
+ add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+ option_category("Restart options");
+ add_option("restart",&restarting,false,"Restart from prior output time.");
+ add_option("restart_time",&initial_time,0.0,"Time to restart from");
+ add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+ option_category("Dumping options");
+ add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+ add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+ option_category("Other options");
+ add_option("perturb",&perturb,"Initial perturbation in velocity");
+ add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+ add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+ add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+ add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+ "Calculate BPE gained from internal energy?");
+ add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+ add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+ add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+ option_category("Filter options");
+ add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+ add_option("f_order",&f_order,2.0,"Filter order");
+ add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+ // 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
+ * can be directly used, but the ones of string-type need further procesing. */
+
+ // adjust temporal values when restarting from dump
+ if (restart_from_dump) {
+ adjust_for_dump(restarting, initial_time, restart_sequence,
+ final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+ }
+
+ // check restart sequence
+ check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+ // parse expansion types
+ parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+ // vector of expansion types
+ grid_type[0] = xgrid_type;
+ grid_type[1] = ygrid_type;
+ grid_type[2] = 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");
+ }
+
+ /* ------------------ Derived parameters --------------------- */
+
+ // Dynamic viscosity
+ mu = visco*rho_0;
+ // Maximum buoyancy frequency (squared) if the initial stratification was stable
+ //N2_max = g*delta_rho/(2*delta_x);
+ // Maximum time step
+ if (dt_max == 0.0) {
+ // if dt_max not given in spins.conf, use this
+ dt_max = 0.1;
+ }
+
+ /* ------------------ Initialize --------------------- */
+
+ // Create an instance of the above class
+ userControl mycode;
+ // Create a flow-evolver that takes its settings from the above class
+ FluidEvolve<userControl> do_stuff(&mycode);
+ // Initialize
+ do_stuff.initialize();
+
+ /* ------------------ Print some parameters --------------------- */
+
+ compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+ double startup_time = compute_start_time - real_start_time;
+
+ if (master()) {
+ fprintf(stdout,"Dam break problem\n");
+ fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+ fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+ fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+ fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+ fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+ fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+ fprintf(stdout,"Max time step: %g\n",dt_max);
+ fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+ }
+
+ /* ------------------ Run --------------------- */
+ // Run until the end of time
+ do_stuff.do_run(final_time);
+ MPI_Finalize();
+ return 0;
+}
diff --git a/src/cases/quadeos_test/spins.conf b/src/cases/quadeos_test/spins.conf
new file mode 100644
index 0000000000000000000000000000000000000000..fe5c278245ccd17011b34d7b7cd6cd4dbdcef1ff
--- /dev/null
+++ b/src/cases/quadeos_test/spins.conf
@@ -0,0 +1,67 @@
+## 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
diff --git a/src/cases/salt_temp/salt_temp.cpp b/src/cases/salt_temp/salt_temp.cpp
index f48746d395f72b0060c1af7c8d125715e39e2f89..99ecd95bf751660bb07d77df3c9ab730f76e8a7b 100644
--- a/src/cases/salt_temp/salt_temp.cpp
+++ b/src/cases/salt_temp/salt_temp.cpp
@@ -244,7 +244,7 @@ class userControl : public BaseCase {
vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
u_f = +rot_f*v;
v_f = -rot_f*u;
- w_f = -g*eqn_of_state(*tracers[TEMP],*tracers[SALT]) / rho_0;
+ w_f = -g*(nleos(*tracers[TEMP],*tracers[SALT]) - rho_0)/rho_0;
*tracers_f[TEMP] = 0;
*tracers_f[SALT] = 0;
}
diff --git a/src/gmres.hpp b/src/gmres.hpp
index 1499ae3234ef06376ffc76b4d12594da2f042e58..863eed30f69ca4096bff355ad44caa9b266b4a5a 100644
--- a/src/gmres.hpp
+++ b/src/gmres.hpp
@@ -269,8 +269,8 @@ template <class Controller> class GMRES_Solver {
/* Perform the inner iteration */
/* First, allocate vectors for residuals and the Krylov basis */
- blitz::Vector<RType> resid_vec(num_its+1);
- blitz::Vector<BType> basis_vec(num_its); // Note, these are base 0
+ blitz::Array<RType,1> resid_vec(num_its+1);
+ blitz::Array<BType,1> basis_vec(num_its); // Note, these are base 0
for (int k = 0; k < num_its; k++) {
resid_vec(k) = ops->alloc_resid();
basis_vec(k) = ops->alloc_basis();
@@ -406,7 +406,7 @@ template <class Controller> class GMRES_Solver {
return my_it-1;
}
- void build_x(BType out_x, blitz::Vector<BType> & basis_vec,
+ void build_x(BType out_x, blitz::Array<BType,1> & basis_vec,
blitz::Array<double,1> & rhs_vec, double starting_norm,
int my_it) {
/* Now, the Krylov basis is in basis_vec, and the proper multiples are in
diff --git a/src/multigrid.cpp b/src/multigrid.cpp
index 7df67096bd518b3e46bf98a718d016a22691b967..46836679ee105f0b58551c511976a9281e5b9433 100644
--- a/src/multigrid.cpp
+++ b/src/multigrid.cpp
@@ -3,7 +3,6 @@
#include <blitz/array.h>
#include <iostream>
-#include <blitz/tinyvec-et.h>
#include "umfpack.h"
#include "timing.hpp"
@@ -1897,14 +1896,14 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
that the problem is really of size (1+size_x*size_z);
this doesn't actually fit in u or f. So, allocate
static blitz-vectors for the data, if necessary. */
- static blitz::Vector<double> extra_u(0), extra_f(0);
+ static blitz::Array<double,1> extra_u(0), extra_f(0);
// fprintf(stderr,"CG: Input residual\n");
// cerr << f;
if (indefinite_problem) {
// fprintf(stderr,"Solving indefinite problem on coarse grid\n");
- if (extra_u.length() != (size_x*size_z+1))
+ if (extra_u.length(firstDim) != (size_x*size_z+1))
extra_u.resize(size_x*size_z+1);
- if (extra_f.length() != (size_x*size_z+1))
+ if (extra_f.length(firstDim) != (size_x*size_z+1))
extra_f.resize(size_x*size_z+1);
// Copy the f-problem to extra_f
diff --git a/systems/belize3.sh b/systems/belize3.sh
new file mode 100644
index 0000000000000000000000000000000000000000..949a94f20210b92d844b8b96a64662c342e8e0c6
--- /dev/null
+++ b/systems/belize3.sh
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+# System-specific settings for belize.math.uwaterloo.ca
+
+CC=gcc
+CXX=g++
+LD=mpic++
+
+# System-specific compiler flags
+SYSTEM_CFLAGS=
+SYSTEM_LDFLAGS=
+SYSTEM_CXXFLAGS=
+
+# Compiler flags for debugging
+DEBUG_CFLAGS="-g -DBZ_DEBUG"
+DEBUG_LDFLAGS=
+
+# Compiler flags for optimization
+OPTIM_CFLAGS="-O3 -ffast-math -march=native"
+OPTIM_LDFLAGS=$OPTIM_CFLAGS
+
+# Compiler flags for extra optimization, such as -ip -ipo on icc
+EXTRA_OPTIM_CFLAGS=
+EXTRA_OPTIM_LDFLAGS=
+
+# Library names/locations/flags for MPI-compilation. This will
+# probably not be necessary on systems with a working mpicc
+# alias
+MPICXX=mpic++
+MPI_CFLAGS=
+MPI_LIB=
+MPI_LIBDIR=
+MPI_INCDIR=
+
+# Library names/locations for LAPACK
+LAPACK_LIB="-llapack -lblas"
+LAPACK_LIBDIR=
+LAPACK_INCDIR=
+
+# Library locations for blitz; leave blank to use system-installed
+# or compiled-with-this-package version
+BLITZ_LIBDIR=
+BLITZ_INCDIR=
+
+# Library locations for fftw
+FFTW_LIBDIR=
+FFTW_INCDIR=
+
+# Library locations for UMFPACK
+UMF_INCDIR=
+UMF_LIBDIR=
+
+# Location/library for BLAS
+BLAS_LIB="-lblas"
+BLAS_LIBDIR=
+BLAS_INCDIR=
+
diff --git a/systems/hood.sh b/systems/hood.sh
index 4cf798bcb78535dd1b8bf54d4b7c5194cbd5ffe2..d6ff86b76cc525171c89bf183f680c089d55f968 100644
--- a/systems/hood.sh
+++ b/systems/hood.sh
@@ -28,9 +28,9 @@ EXTRA_OPTIM_LDFLAGS=$EXTRA_OPTIM_CFLAGS
# alias
MPICXX=icpc
MPI_CFLAGS=
-MPI_LIBDIR="-L${MPI_LIB}"
+MPI_LIBDIR="-L${MPI_LIB?MPI_LIB unset}"
MPI_LIB="-lmpi"
-MPI_INCDIR="-I${MPI_INCLUDE}"
+MPI_INCDIR="-I${MPI_INCLUDE?MPI_INCLUDE unset}"
# Library names/locations for LAPACK
LAPACK_LIB="-lmkl_intel_lp64 -lmkl_core -lmkl_sequential -lpthread"
diff --git a/systems/makemake.sh b/systems/makemake.sh
index b70c4f248a9f34009b78f2c8f8a52326629a8bd3..359f1407c6edfcf7897f659ba16fbdd1cd5ca38d 100755
--- a/systems/makemake.sh
+++ b/systems/makemake.sh
@@ -83,3 +83,18 @@ BLAS_INCDIR:=$BLAS_INCDIR
EOF
+echo "Sanity check - building trivial MPI executable"
+rm -f mpitest.o mpitest.x mpitest.cpp
+cat > mpitest.cpp << EOF
+#include <mpi.h>
+int main(int argc, char ** argv) {
+ MPI_Init(&argc, &argv);
+ MPI_Finalize();
+ return 0;
+}
+EOF
+echo "Executing: $MPICXX $SYSTEM_CFLAGS $SYSTEM_CXXFLAGS $MPI_CFLAGS $MPI_INCDIR -c -o mpitest.o ./mpitest.cpp"
+$MPICXX $SYSTEM_CFLAGS $SYSTEM_CXXFLAGS $MPI_CFLAGS $MPI_INCDIR -c -o mpitest.o ./mpitest.cpp || echo "WARNING: MPI compilation failed!" && echo "MPI compilation success"
+echo "Executing: $LD $SYSTEM_LDFLAGS $MPI_LIBDIR -o mpitest.x mpitest.o $MPI_LIB -lstdc++"
+$LD $SYSTEM_LDFLAGS $MPI_LIBDIR -o mpitest.x mpitest.o $MPI_LIB -lstdc++ || echo "WARNING: MPI linking failed!" && echo "MPI linking success"
+rm -f mpitest.o mpitest.x mpitest.cpp