Commit 31978d73 authored by Christopher Subich's avatar Christopher Subich
Browse files

Merge branch 'master' of git.uwaterloo.ca:ncastrof/SPINS_main

Resolved merge conflicts in Science.hpp, derivatves.cpp, incremented
VERSION by one further patchlevel
parents d63d8368 e4bf9871
......@@ -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)
......
#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);
}
MAJOR_VERSION=2
MINOR_VERSION=1
PATCH_VERSION=7
PATCH_VERSION=8
......@@ -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
......
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment