Commit 14d4ec38 authored by Christopher Subich's avatar Christopher Subich
Browse files

Merge branch 'master' into 'master'

Master

See merge request !6
parents 6fbdee14 ae2e1992
......@@ -110,6 +110,18 @@ void enstrophy_stretch_production(TArrayn::DTArray & enst_prod, TArrayn::DTArray
TArrayn::DTArray & temp2, TArrayn::DTArray & temp3, TArrayn::Grad * gradient_op,
const string * grid_type);
// Q/Second invariant of grad(u,v,w)
void Q_invt(TArrayn::DTArray & Q, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type);
// R/Third invariant of grad(u,v,w)
void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type, bool v_exist);
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
......
#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;
// Q, or the second invariant of grad(u,v,w)
// Defined as u_x * v_y + v_y * w_z + u_x * w_z - u_y * v_x - v_z * w_y - u_z * w_x
// This is used to find coherent structures correlated with eddies in turbulent flow
void Q_invt(TArrayn::DTArray & Q, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
// First term
// Get u_x
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp1,false);
// Get v_y
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// u_x * v_y
Q = temp1 * temp2;
// Second term
// Already have v_y, so get w_z
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// w_z * v_y
Q += temp1 * temp2;
// Third term
// Already have w_z, so get u_z
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// w_z * u_x
Q += temp1 * temp2;
// Fourth term
// Get u_y
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp1,false);
// Get v_x
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// - u_y * v_x
Q -= temp1 * temp2;
// Fifth term
// Get v_z
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get w_y
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// - u_y * v_x
Q -= temp1 * temp2;
// Sixth term
// Get u_z
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get v_x
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// - u_z * w_x
Q -= temp1 * temp2;
}
#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;
// R, the third invariant of grad(u,v,w), which is just det(grad(u,v,w)).
// For 2D defined as:
// R = u_x * w_z - u_z * w_x
// For 3D defined as:
// R = u_x * (v_y * w_z - v_z * w_y) + u_y * (v_z * w_x - v_x * w_z) + u_z * (v_x * w_y - v_y * w_x)
// This is used to find coherent structures correlated with eddies in turbulent flow
void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type, bool v_exist) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
if (!v_exist) {
//2D R
// First term
// Get u_x
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp1,false);
// Get w_z
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// u_x * w_z
R = temp1 * temp2;
// Second term
// Get u_z
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get w_x
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// u_x * w_z
R -= temp1 * temp2;
}
else {
// 3D
// First term
// Get v_y
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp1,false);
// Get w_z
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// v_y * w_z
temp1 = temp1 * temp2;
// Get u_x
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// v_y * w_z * u_x
R = temp1 * temp2;
// Get v_z
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get w_y
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// v_z * w_y
temp1 = temp1 * temp2;
// Get u_x
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// - v_z * w_y * u_x
R -= temp1 * temp2;
// Second term
// Get v_z
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp1,false);
// Get w_x
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// v_z * w_x
temp1 = temp1 * temp2;
// Get u_y
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// + v_z * w_x * u_y
R += temp1 * temp2;
// Get v_x
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp1,false);
// Get w_z
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// v_x * w_z
temp1 = temp1 * temp2;
// Get u_y
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// - v_x * w_z * u_y
R -= temp1 * temp2;
// Third term
// Get v_x
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp1,false);
// Get w_y
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// v_x * w_y
temp1 = temp1 * temp2;
// Get u_z
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// + v_x * w_y * u_z
R += temp1 * temp2;
// Get v_y
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp1,false);
// Get w_x
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// v_y * w_x
temp1 = temp1 * temp2;
// Get u_z
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// - v_y * w_x * u_z
R -= temp1 * temp2;
}
}
MAJOR_VERSION=2
MINOR_VERSION=0
MINOR_VERSION=1
PATCH_VERSION=6
......@@ -46,6 +46,9 @@ bool do_enstrophy; // Do Enstrophy calculation?
bool do_dissipation; // Do Viscous dissipation?
bool do_vort_stretch; // Do vortex stretching/tilting?
bool do_enst_stretch; // Do enstrophy production term from vortex streching/tilting?
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 v_exist; // Does the v field exist?
/* ------------------ Adjust the class --------------------- */
......@@ -78,7 +81,7 @@ class userControl : public BaseCase {
/* Set other things */
double get_visco() const { return visco; }
int get_restart_sequence() const { return plotnum; }
/* Read grid (if mapped) */
bool is_mapped() const { return mapped; }
void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
......@@ -107,13 +110,13 @@ 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 ) {
if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or do_Q_and_R ) {
temp1 = alloc_array(Nx,Ny,Nz);
temp2 = alloc_array(Nx,Ny,Nz);
if ( do_enst_stretch )
if ( do_enst_stretch )
temp3 = alloc_array(Nx,Ny,Nz);
}
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
plotnum = plotnum + step_sequence ) {
......@@ -204,9 +207,11 @@ class userControl : public BaseCase {
}
}
// read in fields (if not already stored in memory)
// read in fields (if not already stored in memory) Nico: Should do_Q go here?
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 ) {
do_enstrophy or do_dissipation or
do_vort_stretch or do_enst_stretch or
do_Q or do_R or do_Q_and_R ) {
// u
init_tracer_restart("u",u);
// v
......@@ -315,6 +320,24 @@ class userControl : public BaseCase {
if (master())
fprintf(stdout,"Completed the write for enst-stretch.%d\n",plotnum);
}
// Calculate Q/second invariant of grad(u,v,w)
if ( do_Q or do_Q_and_R ) {
Q_invt(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max Q: %.6g\n",max_var);
write_array(deriv_var,"Q",plotnum);
if (master()) fprintf(stdout,"Completed the write for Q.%d\n",plotnum);
}
// Calculate R/third invariant of grad(u,v,w)
if ( do_R or do_Q_and_R ) {
R_invt(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type, v_exist);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max R: %.6g\n",max_var);
write_array(deriv_var,"R",plotnum);
if (master()) fprintf(stdout,"Completed the write for R.%d\n",plotnum);
}
}
}
......@@ -380,8 +403,10 @@ int main(int argc, char ** argv) {
add_option("do_dissipation",&do_dissipation,false,"Calculate viscous dissipation?");
add_option("do_vort_stretch",&do_vort_stretch,false,"Calculate vortex stretching?");
add_option("do_enst_stretch",&do_enst_stretch,false,"Calculate enstrophy stretching production?");
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("v_exist",&v_exist,"Does the v field exist?");
// Parse the options from the command line and config file
options_parse(argc,argv);
......
......@@ -34,3 +34,6 @@ do_enstrophy = false
do_dissipation = false
do_vort_stretch = false
do_enst_stretch = false
do_Q = false
do_R = false
do_Q_and_R = true
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