Commit 165b9d88 authored by Nico Castro-Folker's avatar Nico Castro-Folker
Browse files

Update derivatives.cpp to include call for lambda2

parent d9056496
......@@ -32,6 +32,8 @@ const int z_ind = 2;
// physical parameters
double visco; // viscosity (m^2/s)
double rho_0; // reference density (kg/m^3)
double g; // acceleration due to gravity (m/s^2)
// current output number
int plotnum;
......@@ -42,6 +44,7 @@ int final_sequence; // output number to stop taking derivatives
int step_sequence; // step between outputs to take derivatives
bool deriv_x, deriv_y, deriv_z; // which derivatives
bool do_vor_x, do_vor_y, do_vor_z; // Do vorticity calculations?
bool do_barvor; // Do baroclinic vorticity?
bool do_enstrophy; // Do Enstrophy calculation?
bool do_dissipation; // Do Viscous dissipation?
bool do_vort_stretch; // Do vortex stretching/tilting?
......@@ -49,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 --------------------- */
......@@ -59,6 +63,8 @@ class userControl : public BaseCase {
Grad * gradient_op; // gradient operator
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; }
......@@ -110,13 +116,25 @@ 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 )
temp3 = alloc_array(Nx,Ny,Nz);
}
if ( do_barvor ) {
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 ) {
......@@ -206,12 +224,38 @@ class userControl : public BaseCase {
}
}
}
// Baroclinic vorticity
if ( do_barvor ) {
// Store Temperature in T, it is free
init_tracer_restart("t",u);
compute_baroclinic_vort(deriv_var, *temp4, u, gradient_op, grid_type, v_exist);
deriv_var=deriv_var*g;
write_array(deriv_var,"bar",plotnum);
if (master()) fprintf(stdout,"Completed the write for bar.%d\n",plotnum);
if ( v_exist ) {
compute_baroclinic_vort_x(deriv_var, u, gradient_op, grid_type);
deriv_var=deriv_var*g;
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);
if (master()) fprintf(stdout,"Completed the write for bary.%d\n",plotnum);
// Restart u
u = 0;
}
// read in fields (if not already stored in memory) Nico: Should do_Q go here?
// 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 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
......@@ -255,7 +299,7 @@ class userControl : public BaseCase {
write_array(deriv_var,"vortz",plotnum);
if (master())
fprintf(stdout,"Completed the write for vortz.%d\n",plotnum);
}
}
// Calculate Enstrophy
if ( do_enstrophy ) {
......@@ -338,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);
}
}
}
......@@ -385,8 +440,10 @@ 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");
option_category("Derivative options");
add_option("deriv_files",&deriv_filenames,"Derivative filename");
......@@ -399,6 +456,7 @@ int main(int argc, char ** argv) {
add_option("do_vor_x",&do_vor_x,false,"Do the X-component of vorticity?");
add_option("do_vor_y",&do_vor_y,false,"Do the Y-component of vorticity?");
add_option("do_vor_z",&do_vor_z,false,"Do the Z-component of vorticity?");
add_option("do_barvor",&do_barvor,false,"Do the baroclinic 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?");
......@@ -406,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);
......@@ -428,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
......
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