diff --git a/src/Science.hpp b/src/Science.hpp
index 569b9a6c6c494ddbb583efef4e924702b6e43237..4110a58d0829a8b2d315efbdda1ae0fa66156d8b 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -130,9 +130,126 @@ 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)
+*/
+
+
// 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?)
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/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/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/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;
}