Commit 01972559 authored by Donovan Allum's avatar Donovan Allum

Added the baroclinic vorticity as a new output for the derivatives casefile

parent 6fbdee14
......@@ -35,6 +35,14 @@ void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArra
TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type);
// Baroclinic Vorticity
void compute_baroclinic_vort_x(TArrayn::DTArray & barovortx, TArrayn::DTArray & T,
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_baroclinic_vort_y(TArrayn::DTArray & barovorty, TArrayn::DTArray & T,
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp,
TArrayn::DTArray & T, TArrayn::Grad * gradient_op, const string * grid_type, bool v_exist);
// Enstrophy density
void enstrophy_density(TArrayn::DTArray & enst, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
......@@ -134,6 +142,49 @@ inline double eqn_of_state(double T, double S){
// Define a Blitz-friendly operator
BZ_DECLARE_FUNCTION2(eqn_of_state)
// Derivative of the equation of state WRT T
inline double eqn_of_state_dT(double T){
// Returbs the derivative of the eqn_of_state function above WRT T
// temperature T (degrees celsius), S assumed to be zero
// Constants are from table 4 of the above paper, at pressure 0
// (This is appropriate since this is currently an incompressible
// model)
// Kept the same names for these constants as in eqn_of_state
const double c2 = 5.10768e-2; // T term, constant after derivative
const double c4 = -7.40849e-3; // T^2 term, 2T after derivative
const double c5 = -3.01036e-3; // ST term, S after derivative
const double c6 = 3.32267e-5; // T^3 term, 3t^2 after derivative
const double c7 = 3.21931e-5; // ST^2 term, 2ST after derivative
return c2 + 2*c4*T + 3*c6*T*T;
}
// Define a Blitz-friendly operator
BZ_DECLARE_FUNCTION(eqn_of_state_dT)
// Derivative of the equation of state WRT S
inline double eqn_of_state_dS(double T, double S){
// Returns the density anomaly (kg/m^3) for water at the given
// temperature T (degrees celsius) and salinity S (PSU?)
// Constants are from table 4 of the above paper, at pressure 0
// (This is appropriate since this is currently an incompressible
// model)
const double c3 = 8.05999e-1; // S term, constant after derivative
const double c5 = -3.01036e-3; // ST term, T after derivative
const double c7 = 3.21931e-5; // ST^2 term, T^2 after derivative
return c3 + c5*T + c7*T*T;
}
// Define a Blitz-friendly operator
BZ_DECLARE_FUNCTION2(eqn_of_state_dS)
inline double eqn_of_state_t(double T){
// Specialize for freshwater with S=0
return eqn_of_state(T,0.0);
......
#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;
// X-component of the baroclinic vorticity -g*d\rho/dy
void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp,
TArrayn::DTArray & T, TArrayn::Grad * gradient_op, const string * grid_type, bool v_exist ) {
if ( v_exist ) {
compute_baroclinic_vort_x(temp, T, gradient_op, grid_type); }
else {
barovort = 0;
}
barovort = temp*temp;
compute_baroclinic_vort_y(temp, T, gradient_op, grid_type);
barovort = barovort + temp*temp;
barovort = sqrt(barovort);
}
#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;
// X-component of the baroclinic vorticity -g*d\rho/dy
void compute_baroclinic_vort_x(TArrayn::DTArray & barovortx, TArrayn::DTArray & T,
TArrayn::Grad * gradient_op, const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
// Setup for dT/dy
find_expansion(grid_type, expan, "t");
gradient_op->setup_array(&T,expan[0],expan[1],expan[2]);
// get dT/dy
gradient_op->get_dy(&barovortx,false);
// use chainrule to calculate barovortx
barovortx = (-1)*barovortx*eqn_of_state_dT(T);
}
#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;
// Y-component of the baroclinic vorticity g*d\rho/dx=g*d\rho/dT*dT/dx
void compute_baroclinic_vort_y(TArrayn::DTArray & barovorty, TArrayn::DTArray & T,
TArrayn::Grad * gradient_op, const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
// Setup for dT/dx
find_expansion(grid_type, expan, "t");
gradient_op->setup_array(&T,expan[0],expan[1],expan[2]);
// get dT/dx
gradient_op->get_dx(&barovorty,false);
// use chainrule to calculate barovortx
// eqn_of_state_dT(T,0.0) is d\rho/dT
barovorty = barovorty*eqn_of_state_dT(T);
}
......@@ -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?
......@@ -56,6 +59,7 @@ 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
/* Size of domain */
double length_x() const { return Lx; }
......@@ -113,6 +117,10 @@ class userControl : public BaseCase {
if ( do_enst_stretch )
temp3 = alloc_array(Nx,Ny,Nz);
}
if ( do_barvor ) {
temp4 = alloc_array(Nx,Ny,Nz);
}
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
......@@ -203,9 +211,36 @@ 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)
if ( do_vor_x or do_vor_y or do_vor_z or
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);
......@@ -250,7 +285,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 ) {
......@@ -364,6 +399,8 @@ int main(int argc, char ** argv) {
option_category("Physical parameters");
add_option("visco",&visco,0.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");
......@@ -376,6 +413,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?");
......
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