diff --git a/src/Science.hpp b/src/Science.hpp
index 4110a58d0829a8b2d315efbdda1ae0fa66156d8b..24b459c0b039755d14cb37db5eec7fa9ead1b304 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -183,8 +183,6 @@ inline double nleos_inline(double T, double S){
}
BZ_DECLARE_FUNCTION2(nleos_inline)
-
-
void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
TArrayn::DTArray & S);
@@ -229,8 +227,6 @@ 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
@@ -246,6 +242,13 @@ inline double fresh_quad(double T){
BZ_DECLARE_FUNCTION(fresh_quad)
*/
+// lambda2, second eigenvalue of S^2+Omega^2
+void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
+ TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
+ TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
+ const string * grid_type, TArrayn::DTArray & A11, TArrayn::DTArray & A12,
+ TArrayn::DTArray & A13, TArrayn::DTArray & A22, TArrayn::DTArray & A23,
+ TArrayn::DTArray & A33);
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
diff --git a/src/Science/compute_lambda2.cpp b/src/Science/compute_lambda2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abd08d62d940c41c535bfe9a72b3e48cc179ff42
--- /dev/null
+++ b/src/Science/compute_lambda2.cpp
@@ -0,0 +1,161 @@
+#include "../Science.hpp"
+#include <blitz/array.h>
+#include <string>
+#include <math.h>
+
+using namespace Transformer;
+using namespace blitz;
+
+// Calculates the second largest eigenvalue of S^2+Omega^2
+// S and Omega are the symmetric and antisymmetric components of u_i,j, respectively
+// Note: S^2+Omega^2 = 0.5*(u_{i,k}u_{k,j}+u_{k,i}u_{j,k})
+// This is used to find coherent structures correlated with eddies in turbulent flow
+void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
+ TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
+ TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
+ const string * grid_type, TArrayn::DTArray & A11, TArrayn::DTArray & A12,
+ TArrayn::DTArray & A13, TArrayn::DTArray & A22, TArrayn::DTArray & A23,
+ TArrayn::DTArray & A33) {
+ // Set-up
+ S_EXP expan[3];
+ assert(gradient_op);
+ TArrayn::DTArray * A_ref; //This will hold references to the components
+ TArrayn::DTArray * vel_mm; //This will hold the velocities
+ TArrayn::DTArray * vel_nn; //This will hold the velocities
+ std::string vel_labs[3] = {"u","v","w"};
+
+ int ctr = 0;
+ // Construct S^2+Omega^2
+ for (int mm=0; mm<3; mm++){
+ for (int nn=mm; nn<3; nn++){
+ // Choose which component will be constructed
+ if (ctr==0){A_ref = &A11;}
+ else if (ctr==1){A_ref = &A12;}
+ else if (ctr==2){A_ref = &A13;}
+ else if (ctr==3){A_ref = &A22;}
+ else if (ctr==4){A_ref = &A23;}
+ else if (ctr==5){A_ref = &A33;}
+
+ // A switch for the velocities
+ // Useful for computing the u_{i,k}u_{j,k} terms
+
+ if (mm==0){vel_mm = &u;}
+ else if (mm==1){vel_mm = &v;}
+ else if(mm==2){vel_mm = &w;}
+ if (nn==0){vel_nn = &u;}
+ else if (nn==1){vel_nn = &v;}
+ else if (nn==2){vel_nn = &w;}
+
+ // u_{i,k}u_{k,j}
+ // Get u_{i,x}
+ find_expansion(grid_type, expan, vel_labs[mm]);
+ gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+ gradient_op->get_dx(&temp1,false);
+ // Get u_{x,j}
+ find_expansion(grid_type, expan, "u");
+ gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+ if (nn==0) { gradient_op->get_dx(&temp2,false); }
+ else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+ else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+ *A_ref = temp1 * temp2;
+ // Get u_{i,y}
+ find_expansion(grid_type, expan, vel_labs[mm]);
+ gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+ gradient_op->get_dy(&temp1,false);
+ // Get u_{y,j}
+ find_expansion(grid_type, expan, "v");
+ gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+ if (nn==0) { gradient_op->get_dx(&temp2,false); }
+ else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+ else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+ *A_ref = temp1 * temp2;
+ // Get u_{i,z}
+ find_expansion(grid_type, expan, vel_labs[mm]);
+ gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+ gradient_op->get_dz(&temp1,false);
+ // Get u_{z,j}
+ find_expansion(grid_type, expan, "w");
+ gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+ if (nn==0) { gradient_op->get_dx(&temp2,false); }
+ else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+ else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+ *A_ref = temp1 * temp2;
+
+
+ // u_{k,i}u_{j,k}
+ // Get u_{x,i}
+ find_expansion(grid_type, expan, "u");
+ gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+ if (mm==0) { gradient_op->get_dx(&temp1,false); }
+ else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+ else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+ // Get u_{j,x}
+ find_expansion(grid_type, expan, vel_labs[nn]);
+ gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+ gradient_op->get_dx(&temp2,false);
+ *A_ref = temp1 * temp2;
+ // Get u_{y,i}
+ find_expansion(grid_type, expan, "v");
+ gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+ if (mm==0) { gradient_op->get_dx(&temp1,false); }
+ else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+ else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+ // Get u_{j,y}
+ find_expansion(grid_type, expan, vel_labs[nn]);
+ gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+ gradient_op->get_dy(&temp2,false);
+ *A_ref = temp1 * temp2;
+ // Get u_{z,i}
+ find_expansion(grid_type, expan, "w");
+ gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+ if (mm==0) { gradient_op->get_dx(&temp1,false); }
+ else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+ else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+ // Get u_{j,z}
+ find_expansion(grid_type, expan, vel_labs[nn]);
+ gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+ gradient_op->get_dz(&temp2,false);
+ *A_ref = temp1 * temp2;
+
+ *A_ref = (*A_ref)*0.5;
+ ctr++;
+ }
+ }
+
+ // Now for the eigenvalue algorithm from:
+ // https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
+
+ // In the above article, we transform A into B=(1/p)*(A-qI)
+ // Then, lambda2 = q + 2*p*cos(arccos(det(B)/2)/3 + k*pi/3), k=0,1,2
+ // Define q
+ temp1 = (A11+A22+A33)/3.0;
+ // Define p
+ temp2 = pow2(A11-temp1) + pow2(A22-temp1) + pow2(A33-temp1) +2.0*(pow2(A12)+pow2(A13)+pow2(A23));
+ temp2 = sqrt(temp2/6.0);
+
+ // Transform A
+ A11 = (A11-temp1)/temp2;
+ A22 = (A22-temp1)/temp2;
+ A33 = (A33-temp1)/temp2;
+ A12 = A12/temp2;
+ A13 = A13/temp2;
+ A23 = A23/temp2;
+
+ // Calculate the determinant, using lambda2 as a dummy array
+ lambda2 = 0.5*(A11*(A22*A33-A23*A23)-A12*(A12*A33-A13*A23)+A13*(A12*A23-A13*A22));
+
+ // Since det(B)/2 feeds into acos, make sure it's between -1 and 1
+ for (int i = lambda2.lbound(blitz::firstDim); i <= lambda2.ubound(blitz::firstDim); i++) { // outer loop over slowest-varying dimension
+ for (int k = lambda2.lbound(blitz::thirdDim); k <= lambda2.ubound(blitz::thirdDim); k++) { // middle loop over next-slowest varying dimension
+ for (int j = lambda2.lbound(blitz::secondDim); j <= lambda2.ubound(blitz::secondDim); j++) { // inner loop over fastest varying dimeion
+ if (isnan(lambda2(i,j,k))) {lambda2(i,j,k) = 0;}
+ else if (lambda2(i,j,k) < -1.0) {lambda2(i,j,k) = -1.0;}
+ else if (lambda2(i,j,k) > 1.0) {lambda2(i,j,k) = 1.0;}
+ // first conditional works as a safety check because the det is NaN iff temp2=0; in this case the 2nd eig is temp1
+ }}}
+
+ // Now calculate the actual lamba2.
+ lambda2 = temp1+2.0*temp2*cos(acos(lambda2)/3.0+4.0*M_PI/3.0);
+
+
+}
diff --git a/src/VERSION b/src/VERSION
index 882b025b720b8ddb68732c2992c7484ee47f79f2..ae8b19475b67f0169e4f6057b5c5728ce674b540 100644
--- a/src/VERSION
+++ b/src/VERSION
@@ -1,3 +1,3 @@
MAJOR_VERSION=2
MINOR_VERSION=1
-PATCH_VERSION=7
+PATCH_VERSION=8
diff --git a/src/cases/derivatives/derivatives.cpp b/src/cases/derivatives/derivatives.cpp
index 21657117390cfa421e678a6415cd0468558fe3d6..e34cf225c1f576bd8c7423d00db8723d79d3b4fb 100644
--- a/src/cases/derivatives/derivatives.cpp
+++ b/src/cases/derivatives/derivatives.cpp
@@ -52,6 +52,7 @@ bool do_enst_stretch; // Do enstrophy production term from vortex
bool do_Q; // Do second invariant of grad(u,v,w)?
bool do_R; // Do third invariant/det of grad(u,v,w)?
bool do_Q_and_R; // Do the second and third invariants?
+bool do_lambda2; // Do lambda2/Hussain's Lambda?
bool v_exist; // Does the v field exist?
/* ------------------ Adjust the class --------------------- */
@@ -63,6 +64,7 @@ class userControl : public BaseCase {
DTArray deriv_var; // array for derivative
DTArray *temp1, *temp2, *temp3; // arrays for vortex stretching / enstrophy production
DTArray *temp4; // array for storing temporary array for baroclinic vorticity
+ DTArray *A11, *A12, *A13, *A22, *A23, *A33; //Arrays for components of S^2+Omega^2
/* Size of domain */
double length_x() const { return Lx; }
@@ -114,7 +116,7 @@ class userControl : public BaseCase {
//string prev_deriv, base_field;
vector<string> fields; // vector of fields to take derivatives
split(deriv_filenames.c_str(), ' ', fields); // populate that vector
- if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or do_Q_and_R ) {
+ if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or do_Q_and_R or do_lambda2 ) {
temp1 = alloc_array(Nx,Ny,Nz);
temp2 = alloc_array(Nx,Ny,Nz);
if ( do_enst_stretch )
@@ -125,6 +127,14 @@ class userControl : public BaseCase {
temp4 = alloc_array(Nx,Ny,Nz);
}
+ if ( do_lambda2 ) {
+ A11 = alloc_array(Nx,Ny,Nz);
+ A12 = alloc_array(Nx,Ny,Nz);
+ A13 = alloc_array(Nx,Ny,Nz);
+ A22 = alloc_array(Nx,Ny,Nz);
+ A23 = alloc_array(Nx,Ny,Nz);
+ A33 = alloc_array(Nx,Ny,Nz);
+ }
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
plotnum = plotnum + step_sequence ) {
@@ -231,7 +241,7 @@ class userControl : public BaseCase {
write_array(deriv_var,"barx",plotnum);
if (master()) fprintf(stdout,"Completed the write for barx.%d\n",plotnum);
}
-
+
compute_baroclinic_vort_y(deriv_var, u, gradient_op, grid_type);
deriv_var=deriv_var*g;
write_array(deriv_var,"bary",plotnum);
@@ -245,7 +255,7 @@ class userControl : public BaseCase {
if ( do_vor_x or do_vor_y or do_vor_z or
do_enstrophy or do_dissipation or
do_vort_stretch or do_enst_stretch or
- do_Q or do_R or do_Q_and_R ) {
+ do_Q or do_R or do_Q_and_R or do_lambda2 ) {
// u
init_tracer_restart("u",u);
// v
@@ -372,6 +382,17 @@ class userControl : public BaseCase {
write_array(deriv_var,"R",plotnum);
if (master()) fprintf(stdout,"Completed the write for R.%d\n",plotnum);
}
+
+ // Calculate lambda2/second eigenvalue of S^2+Omega^2
+ if ( do_lambda2 ){
+ compute_lambda2( deriv_var, u, v, w, *temp1, *temp2, gradient_op,
+ grid_type, *A11, *A12, *A13, *A22, *A23, *A33);
+ double max_var = psmax(max(abs(deriv_var)));
+ if (master()) fprintf(stdout,"Max lam2: %.6g\n",max_var);
+ write_array(deriv_var,"lam2",plotnum);
+ if (master()) fprintf(stdout,"Completed the write for lam2.%d\n",plotnum);
+
+ }
}
}
@@ -419,7 +440,7 @@ int main(int argc, char ** argv) {
add_option("type_z",&zgrid_type,"Grid type in Z");
option_category("Physical parameters");
- add_option("visco",&visco,0.0,"Viscosity");
+ add_option("visco",&visco,1.0,"Viscosity");
add_option("rho_0",&rho_0,1000.0,"Reference Density");
add_option("g",&g,9.81,"Acceleration due to gravity");
@@ -443,6 +464,7 @@ int main(int argc, char ** argv) {
add_option("do_Q",&do_Q,false,"Calculate Q?");
add_option("do_R",&do_R,false,"Calculate R?");
add_option("do_Q_and_R",&do_Q_and_R,false,"Calculate Q and R?");
+ add_option("do_lambda2",&do_lambda2,false,"Calculate Lambda2?");
add_option("v_exist",&v_exist,"Does the v field exist?");
// Parse the options from the command line and config file
options_parse(argc,argv);
@@ -465,6 +487,17 @@ int main(int argc, char ** argv) {
fprintf(stdout,"Simulation is 2 dimensional, "
"Ly has been changed to 1.0 for normalization.\n");
}
+ if (visco==1.0){
+ if (master())
+ fprintf(stdout,"You may have forgotten to specify viscosity, "
+ "Using default value visco = 1.\n");
+ }
+
+ if (rho_0==1000.0){
+ if (master())
+ fprintf(stdout,"You may have forgotten to specify reference density, "
+ "Using default value rho_0 = 1000.\n");
+ }
/* ------------------ Do stuff --------------------- */
userControl mycode; // Create an instantiated object of the above class
diff --git a/src/cases/derivatives/spins.conf b/src/cases/derivatives/spins.conf
index 352e07a7d058b22f19f66eb1deb0da4c008cb935..8b53d2d202b4df0f4f8e2952e8768967c15ba688 100644
--- a/src/cases/derivatives/spins.conf
+++ b/src/cases/derivatives/spins.conf
@@ -40,4 +40,5 @@ do_enst_stretch = false
do_barvor = false
do_Q = false
do_R = false
-do_Q_and_R = true
+do_Q_and_R = false
+do_lambda2 = false