Commit a63ba535 authored by David Deepwell's avatar David Deepwell

Add vortex stretching and enstrophy stretching production terms

Vortex stretching is one of the dominant terms in both the
vorticity equation and enstrophy equation. They also help in
understanding the vortex dynamics and how energy at large scales
cascades into small scales.

The three components of vortex stretching (omega dot grad) u
and the enstrophy stretching production (omega_i omega_j S_ij)
are now included and can be computed with the derivatives
case file. Files for vortex stretching terms are called
vort-stretch<dim> where <dim> is one of x,y,z. The enstrophy
stretching production file is called enst-stretch.

As these quantities are the product of velocity gradients, care
must be made to ensure adequate resolution is present.

Update version to 2.0.3
parent 220ef6e5
......@@ -93,6 +93,23 @@ void bottom_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & Hprime,
TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
const double visco);
// Vortex stretching/tilting
void vortex_stretch_x(TArrayn::DTArray & vort_stretch, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op, const string * grid_type);
void vortex_stretch_y(TArrayn::DTArray & vort_stretch, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op, const string * grid_type);
void vortex_stretch_z(TArrayn::DTArray & vort_stretch, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op, const string * grid_type);
// Enstrophy production via vortex stretching/tilting
void enstrophy_stretch_production(TArrayn::DTArray & enst_prod, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::DTArray & temp3, TArrayn::Grad * gradient_op,
const string * grid_type);
// 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;
// Enstrophy production from vortex stetching/tilting
// Defined as omega_i * omega_j * S_{ij}
// or as omega_i * omega_j * du_i/dx_j
// where index notation has been used (See Davidson)
// This is the vorticity dotted with the vorticity stretching term in the vorticity equation
void enstrophy_stretch_production(TArrayn::DTArray & enst_prod, TArrayn::DTArray & u,
TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
TArrayn::DTArray & temp2, TArrayn::DTArray & temp3, TArrayn::Grad * gradient_op,
const string * grid_type) {
// Set-up
assert(gradient_op);
// x-dimension contribution
vortex_stretch_x(temp3, u, v, w, temp1, temp2, gradient_op, grid_type);
compute_vort_x(temp1, v, w, gradient_op, grid_type);
// omega_x * vortex_stretch_x
enst_prod = temp1 * temp3;
// y-dimension contribution
vortex_stretch_y(temp3, u, v, w, temp1, temp2, gradient_op, grid_type);
compute_vort_y(temp1, u, w, gradient_op, grid_type);
// omega_y * vortex_stretch_y
enst_prod += temp1 * temp3;
// z-dimension contribution
vortex_stretch_z(temp3, u, v, w, temp1, temp2, gradient_op, grid_type);
compute_vort_z(temp1, u, v, gradient_op, grid_type);
// omega_z * vortex_stretch_z
enst_prod += temp1 * temp3;
}
#include "../Science.hpp"
using namespace Transformer;
// x component of vortex stretching (omega dot grad) u
void vortex_stretch_x(TArrayn::DTArray & vort_stretch, 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);
// x-vorticity
compute_vort_x(temp1, v, w, gradient_op, grid_type);
// get du/dx
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// omega_x * du/dx
vort_stretch = temp1 * temp2;
// y-vorticity
compute_vort_y(temp1, u, w, gradient_op, grid_type);
// du/dy
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// omega_y * du/dy
vort_stretch += temp1 * temp2;
// z-vorticity
compute_vort_z(temp1, u, v, gradient_op, grid_type);
// du/dz
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// omega_z * du/dz
vort_stretch += temp1 * temp2;
}
#include "../Science.hpp"
using namespace Transformer;
// y component of vortex stretching (omega dot grad) v
void vortex_stretch_y(TArrayn::DTArray & vort_stretch, 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);
// x-vorticity
compute_vort_x(temp1, v, w, gradient_op, grid_type);
// dv/dx
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// omega_x * dv/dx
vort_stretch = temp1 * temp2;
// y-vorticity
compute_vort_y(temp1, u, w, gradient_op, grid_type);
// dv/dy
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// omega_y * dv/dy
vort_stretch += temp1 * temp2;
// z-vorticity
compute_vort_z(temp1, u, v, gradient_op, grid_type);
// dv/dz
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// omega_z * dv/dz
vort_stretch += temp1 * temp2;
}
#include "../Science.hpp"
using namespace Transformer;
// z component of vortex stretching (omega dot grad) w
void vortex_stretch_z(TArrayn::DTArray & vort_stretch, 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);
// x-vorticity
compute_vort_x(temp1, v, w, gradient_op, grid_type);
// dw/dx
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp2,false);
// omega_x * dw/dx
vort_stretch = temp1 * temp2;
// y-vorticity
compute_vort_y(temp1, u, w, gradient_op, grid_type);
// dw/dy
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dy(&temp2,false);
// omega_y * dw/dy
vort_stretch += temp1 * temp2;
// z-vorticity
compute_vort_z(temp1, u, v, gradient_op, grid_type);
// dw/dz
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp2,false);
// omega_z * dw/dz
vort_stretch += temp1 * temp2;
}
MAJOR_VERSION=2
MINOR_VERSION=0
PATCH_VERSION=2
PATCH_VERSION=3
......@@ -44,6 +44,8 @@ bool deriv_x, deriv_y, deriv_z; // which derivatives
bool do_vor_x, do_vor_y, do_vor_z; // Do vorticity calculations?
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 v_exist; // Does the v field exist?
/* ------------------ Adjust the class --------------------- */
......@@ -53,6 +55,7 @@ class userControl : public BaseCase {
/* Initialize things */
Grad * gradient_op; // gradient operator
DTArray deriv_var; // array for derivative
DTArray *temp1, *temp2, *temp3; // arrays for vortex stretching / enstrophy production
/* Size of domain */
double length_x() const { return Lx; }
......@@ -104,6 +107,12 @@ 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 ) {
temp1 = alloc_array(Nx,Ny,Nz);
temp2 = alloc_array(Nx,Ny,Nz);
if ( do_enst_stretch )
temp3 = alloc_array(Nx,Ny,Nz);
}
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
......@@ -196,7 +205,8 @@ class userControl : public BaseCase {
}
// read in fields (if not already stored in memory)
if ( do_vor_x or do_vor_y or do_vor_z or do_enstrophy or do_dissipation) {
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 ) {
// u
init_tracer_restart("u",u);
// v
......@@ -215,7 +225,7 @@ class userControl : public BaseCase {
}
// X-component of vorticity
if (do_vor_x) {
if ( do_vor_x ) {
compute_vort_x(deriv_var, v, w, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max X-vorticity: %.6g\n",max_var);
......@@ -224,7 +234,7 @@ class userControl : public BaseCase {
fprintf(stdout,"Completed the write for vortx.%d\n",plotnum);
}
// Y-component of vorticity
if (do_vor_y) {
if ( do_vor_y ) {
compute_vort_y(deriv_var, u, w, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max Y-vorticity: %.6g\n",max_var);
......@@ -233,7 +243,7 @@ class userControl : public BaseCase {
fprintf(stdout,"Completed the write for vorty.%d\n",plotnum);
}
// Z-component of vorticity
if (do_vor_z) {
if ( do_vor_z ) {
compute_vort_z(deriv_var, u, v, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max Z-vorticity: %.6g\n",max_var);
......@@ -270,6 +280,41 @@ class userControl : public BaseCase {
if (master())
fprintf(stdout,"Completed the write for diss.%d\n",plotnum);
}
// Calculate vortex stretching term
if ( do_vort_stretch ) {
// x component
vortex_stretch_x(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max vortex stretch in x: %.6g\n",max_var);
write_array(deriv_var,"vort-stretchx",plotnum);
if (master())
fprintf(stdout,"Completed the write for vort-stretchx.%d\n",plotnum);
// y component
vortex_stretch_y(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type);
max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max vortex stretch in y: %.6g\n",max_var);
write_array(deriv_var,"vort-stretchy",plotnum);
if (master())
fprintf(stdout,"Completed the write for vort-stretchy.%d\n",plotnum);
// z component
vortex_stretch_z(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type);
max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max vortex stretch in z: %.6g\n",max_var);
write_array(deriv_var,"vort-stretchz",plotnum);
if (master())
fprintf(stdout,"Completed the write for vort-stretchz.%d\n",plotnum);
}
// Calculate enstrophy production via vortex stretching/tilting
if ( do_enst_stretch ) {
enstrophy_stretch_production(deriv_var, u, v, w, *temp1, *temp2, *temp3, gradient_op, grid_type);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max enstrophy stretch production: %.6g\n",max_var);
write_array(deriv_var,"enst-stretch",plotnum);
if (master())
fprintf(stdout,"Completed the write for enst-stretch.%d\n",plotnum);
}
}
}
......@@ -333,6 +378,8 @@ int main(int argc, char ** argv) {
add_option("do_vor_z",&do_vor_z,false,"Do the Z-component of vorticity?");
add_option("do_enstrophy",&do_enstrophy,false,"Calculate enstrophy?");
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("v_exist",&v_exist,"Does the v field exist?");
// Parse the options from the command line and config file
......
......@@ -13,6 +13,7 @@ type_x = FREE_SLIP
type_y = FOURIER
type_z = NO_SLIP
mapped_grid = false
v_exist = false
# Physical Parameters (for dissipation)
visco = 2e-6
......@@ -31,4 +32,5 @@ do_vor_y = true
do_vor_z = false
do_enstrophy = false
do_dissipation = false
v_exist = false
do_vort_stretch = false
do_enst_stretch = false
Markdown is supported
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