Commit 92b79931 authored by David Deepwell's avatar David Deepwell Committed by Christopher Subich
Browse files

Add case for computing derivatives and vorticity

A new case file, derivatives.cpp, gives a means to compute
derivatives of any field or to compute vorticity components.
The vorticity calculation in Science is completely re-written
so-as to also work for mapped grids.

The derivatives file is built to compute secondary variables
after a run has already been completed. This uses spins' built-in
derivative toolkit allowing multiple variables to be calculated
in parallel. A derivative field (one that is the derivative of
another) is denoted by an underscore followed by the direction
the derivative was taken in (ie. u_x is the x derivative of u).
parent 075caad6
......@@ -330,6 +330,14 @@ void BaseCase::write_chain(const char *filename, DTArray & val, int Iout, int Jo
fclose(fid);
}
/* Read grid from regular output */
void BaseCase::init_grid_restart(const std::string & component,
const std::string & filename, DTArray & grid) {
if (master()) fprintf(stdout,"Reading %s from %s\n",component.c_str(),filename.c_str());
read_array(grid,filename.c_str(),size_x(),size_y(),size_z());
return;
}
/* Check and dump */
void BaseCase::check_and_dump(double clock_time, double real_start_time,
......
......@@ -97,6 +97,8 @@ class BaseCase {
virtual void init_vels_dump(DTArray & u, DTArray & v, DTArray & w);
virtual void init_tracer_restart(const std::string & field, DTArray & the_tracer);
virtual void init_tracer_dump(const std::string & field, DTArray & the_tracer);
virtual void init_grid_restart(const std::string & component,
const std::string & filename, DTArray & grid);
virtual void init_tracer(int t_num, DTArray & tracer) {
assert(0 && "init_tracer not implemented");
......
......@@ -154,123 +154,79 @@ void read_2d_slice(Array<double,3> & fillme, const char * filename,
delete sliced;
}
// X-component of vorticity
void compute_vort_x(TArrayn::DTArray & vortx, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
void vorticity(TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w,
TArrayn::DTArray * & vor_x, TArrayn::DTArray * & vor_y,
TArrayn::DTArray * & vor_z, double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
static int Nx = 0, Ny = 0, Nz = 0;
static Transformer::Trans1D trans_x(szx,szy,szz,firstDim,
(DIM_X == PERIODIC ? FOURIER : REAL)),
trans_y(szx,szy,szz,secondDim,
(DIM_Y == PERIODIC ? FOURIER : REAL)),
trans_z (szx,szy,szz,thirdDim,
(DIM_Z == PERIODIC ? FOURIER : REAL));
static blitz::TinyVector<int,3>
local_lbounds(alloc_lbound(szx,szy,szz)),
local_extent(alloc_extent(szx,szy,szz));
static blitz::GeneralArrayStorage<3>
local_storage(alloc_storage(szx,szy,szz));
static DTArray vort_x(local_lbounds,local_extent,local_storage),
vort_y(local_lbounds,local_extent,local_storage),
vort_z(local_lbounds,local_extent,local_storage),
temp_a(local_lbounds,local_extent,local_storage);
/* Initialization */
if (Nx == 0 || Ny == 0 || Nz == 0) {
Nx = szx; Ny = szy; Nz = szz;
}
assert (Nx == szx && Ny == szy && Nz == szz);
/* x-vorticity is w_y - v_z */
vort_x = 0;
if (szy > 1) { // w_y
if (DIM_X == PERIODIC) {
deriv_fft(w,trans_y,temp_a);
vort_x = temp_a*(2*M_PI/Ly);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(w,trans_y,temp_a);
vort_x = temp_a*(M_PI/Ly);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(w,trans_y,temp_a);
vort_x = temp_a*(-2/Ly);
}
}
if (szz > 1) { // v_z
if (DIM_Z == PERIODIC) {
deriv_fft(v,trans_z,temp_a);
vort_x -= temp_a*(2*M_PI/Lz);
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(v,trans_z,temp_a);
vort_x -= temp_a*(M_PI/Lz);
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(v,trans_z,temp_a);
vort_x -= temp_a*(-2/Lz);
}
}
// y-vorticity is u_z - w_x
vort_y = 0;
if (szz > 1) { // u_z
if (DIM_Z == PERIODIC) {
deriv_fft(u,trans_z,temp_a);
vort_y = temp_a*(2*M_PI/Lz);
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(u,trans_z,temp_a);
vort_y = temp_a*(M_PI/Lz);
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(u,trans_z,temp_a);
vort_y = temp_a*(-2/Lz);
}
}
if (szx > 1) { // w_x
if (DIM_X == PERIODIC) {
deriv_fft(w,trans_x,temp_a);
vort_y -= temp_a*(2*M_PI/Lx);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(w,trans_x,temp_a);
vort_y -= temp_a*(M_PI/Lx);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(w,trans_x,temp_a);
vort_y -= temp_a*(-2/Lx);
}
}
// And finally, vort_z is v_x - u_y
vort_z = 0;
if (szx > 1) { // v_x
if (DIM_X == PERIODIC) {
deriv_fft(v,trans_x,temp_a);
vort_z = temp_a*(2*M_PI/Lx);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(v,trans_x,temp_a);
vort_z = temp_a*(M_PI/Lx);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(v,trans_x,temp_a);
vort_z = temp_a*(-2/Lx);
}
}
if (szy > 1) { // u_y
if (DIM_Y == PERIODIC) {
deriv_fft(u,trans_y,temp_a);
vort_z -= temp_a*(2*M_PI/Ly);
} else if (DIM_Y == FREE_SLIP) {
deriv_dct(u,trans_y,temp_a);
vort_z -= temp_a*(M_PI/Ly);
} else {
assert(DIM_Y == NO_SLIP);
deriv_cheb(u,trans_y,temp_a);
vort_z -= temp_a*(-2/Ly);
}
}
vor_x = &vort_x;
vor_y = &vort_y;
vor_z = &vort_z;
return;
// Setup for dv/dz
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
// get dv/dz
gradient_op->get_dz(&vortx,false);
// Invert to get the negative
vortx = (-1)*vortx;
// Setup for dw/dy
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
// get dw/dy, and add to vortx
gradient_op->get_dy(&vortx,true);
}
// Y-component of vorticity
void compute_vort_y(TArrayn::DTArray & vorty, TArrayn::DTArray & u, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
// Setup for dw/dx
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
// get dw/dx
gradient_op->get_dx(&vorty,false);
// Invert to get the negative
vorty = (-1)*vorty;
// Setup for du/dz
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
// get du/dz, and add to vorty
gradient_op->get_dz(&vorty,true);
}
// Z-component of vorticity
void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::Grad * gradient_op, const string * grid_type) {
// Set-up
S_EXP expan[3];
assert(gradient_op);
// Setup for du/dy
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
// get du/dy
gradient_op->get_dy(&vortz,false);
// Invert to get the negative
vortz = (-1)*vortz;
// Setup for dv/dx
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
// get dv/dx, and add to vortz
gradient_op->get_dx(&vortz,true);
}
void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArrayn::DTArray & vortz,
TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type) {
// compute each component
compute_vort_x(vortx, v, w, gradient_op, grid_type);
compute_vort_y(vorty, u, w, gradient_op, grid_type);
compute_vort_z(vortz, u, v, gradient_op, grid_type);
}
// Global arrays to store quadrature weights
......@@ -361,3 +317,72 @@ const blitz::Array<double,1> * get_quad_z() {
}
return &_quadw_z;
}
// function to parse the expansion types
void find_expansion(const string * grid_type, S_EXP * expan, string deriv_filename) {
const int x_ind = 0;
const int y_ind = 1;
const int z_ind = 2;
string prev_deriv, base_field;
// check if field is a derivative field
bool input_deriv = false; // assume it's not a derivative
int var_len = deriv_filename.length();
if ( var_len > 2 ) {
if ( deriv_filename.substr(var_len-2,1) == "_" ) {
// if second last char is an underscore then its a derivative field
input_deriv = true;
prev_deriv = deriv_filename.substr(var_len-1,1); // the completed derivative
base_field = deriv_filename.substr(0,var_len-2); // the differentiated field
}
}
for ( int nn = 0; nn <= 2; nn++ ) {
if (grid_type[nn] == "FOURIER") { expan[nn] = FOURIER; }
else if (grid_type[nn] == "NO_SLIP") { expan[nn] = CHEBY; }
else if (grid_type[nn] == "FREE_SLIP") {
// setup for a first derivative
if ( deriv_filename == "u" or base_field == "u" ) {
if ( nn == x_ind ) { expan[nn] = SINE; }
else if ( nn == y_ind ) { expan[nn] = COSINE; }
else if ( nn == z_ind ) { expan[nn] = COSINE; }
}
else if ( deriv_filename == "v" or base_field == "v" ) {
if ( nn == x_ind ) { expan[nn] = COSINE; }
else if ( nn == y_ind ) { expan[nn] = SINE; }
else if ( nn == z_ind ) { expan[nn] = COSINE; }
}
else if ( deriv_filename == "w" or base_field == "w") {
if ( nn == x_ind ) { expan[nn] = COSINE; }
else if ( nn == y_ind ) { expan[nn] = COSINE; }
else if ( nn == z_ind ) { expan[nn] = SINE; }
}
else {
if ( nn == x_ind ) { expan[nn] = COSINE; }
else if ( nn == y_ind ) { expan[nn] = COSINE; }
else if ( nn == z_ind ) { expan[nn] = COSINE; }
}
}
}
// adjust if input field is a derivative field
if ( input_deriv == true ) {
if ( prev_deriv == "x" ) { expan[x_ind] = swap_trig(expan[x_ind]); }
else if ( prev_deriv == "y" ) { expan[y_ind] = swap_trig(expan[y_ind]); }
else if ( prev_deriv == "z" ) { expan[z_ind] = swap_trig(expan[z_ind]); }
}
}
// function to switch trig functions
S_EXP swap_trig( S_EXP the_exp ) {
if ( the_exp == SINE ) {
return COSINE; }
else if ( the_exp == COSINE ) {
return SINE; }
else if ( the_exp == FOURIER ) {
return FOURIER; }
else if ( the_exp == CHEBY ) {
return CHEBY; }
else {
MPI_Finalize(); exit(1); // stop
}
}
......@@ -24,13 +24,16 @@ void read_2d_slice(blitz::Array<double,3> & fillme, const char * filename,
void read_2d_restart(blitz::Array<double,3>& fillme, const char* filename,
int Nx, int Ny);
/* Compute vorticity */
void vorticity(TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::DTArray * & w_x, TArrayn::DTArray * & w_y,
TArrayn::DTArray * & w_z, double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
// Vorticity
void compute_vort_x(TArrayn::DTArray & vortx, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_vort_y(TArrayn::DTArray & vorty, TArrayn::DTArray & u, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::Grad * gradient_op, const string * grid_type);
void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArrayn::DTArray & vortz,
TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::Grad * gradient_op, const string * grid_type);
// Quadrature weights
......@@ -42,6 +45,11 @@ const blitz::Array<double,1> * get_quad_x();
const blitz::Array<double,1> * get_quad_y();
const blitz::Array<double,1> * get_quad_z();
// find which expansion to use based on field and boundary conditions
void find_expansion(const string * grid_type, Transformer::S_EXP * expan, string deriv_filename);
// switch trig function
Transformer::S_EXP swap_trig( Transformer::S_EXP the_exp );
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
......
/* Derivative case file for computing derivatives of existing fields */
// only able to return first and second derivatives
/* ------------------ Top matter --------------------- */
// Required headers
#include "../../BaseCase.hpp" // contains default class
#include "../../Options.hpp" // config-file parser
#include "../../Science.hpp" // Science content
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::thirdIndex kk;
/* ------------------ Define parameters --------------------- */
// Grid scales
double Lx, Ly, Lz; // Grid lengths (m)
int Nx, Ny, Nz; // Number of points in x, y, z
// Mapped grid?
bool mapped;
// Grid types
DIMTYPE intype_x, intype_y, intype_z;
string grid_type[3];
// Expansion types and indices
S_EXP expan[3];
const int x_ind = 0;
const int y_ind = 1;
const int z_ind = 2;
// physical parameters
double visco; // viscosity (m^2/s)
double rho_0; // reference density (kg/m^3)
// current output number
int plotnum;
// Derivative options
string deriv_filenames; // file name to take derivative of
int start_sequence; // output number to start taking derivatives at
int final_sequence; // output number to stop taking derivatives at
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_enstrophy; // Do Enstrophy calculation?
bool do_dissipation; // Do Viscous dissipation?
bool v_exist; // Does the v field exist?
/* ------------------ Adjust the class --------------------- */
class userControl : public BaseCase {
public:
/* Initialize things */
Grad * gradient_op; // gradient operator
DTArray deriv_var; // array for derivative
/* Size of domain */
double length_x() const { return Lx; }
double length_y() const { return Ly; }
double length_z() const { return Lz; }
/* Resolution in X, Y, and Z */
int size_x() const { return Nx; }
int size_y() const { return Ny; }
int size_z() const { return Nz; }
/* Set expansions (FREE_SLIP, NO_SLIP (in vertical) or PERIODIC) */
DIMTYPE type_x() const { return intype_x; }
DIMTYPE type_y() const { return intype_y; }
DIMTYPE type_z() const { return intype_z; }
/* Record the gradient-taking object */
void set_grad(Grad * in_grad) { gradient_op = in_grad; }
/* 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) {
init_grid_restart("x","xgrid",xg);
if ( Ny > 1 )
init_grid_restart("y","ygrid",yg);
init_grid_restart("z","zgrid",zg);
}
/* function for splitting a string - to parse deriv_filenames */
void split(const string &s, char delim, vector<string> &elems) {
stringstream ss(s);
string item;
while (getline(ss, item, delim)) {
elems.push_back(item);
}
}
/* Read fields and do derivatives */
void init_vels(DTArray & u, DTArray & v, DTArray & w) {
// set-up
assert(gradient_op);
char filename[100];
bool saved_v = false;
bool saved_w = false;
//string prev_deriv, base_field;
vector<string> fields; // vector of fields to take derivatives
split(deriv_filenames.c_str(), ' ', fields); // populate that vector
// Compute derivatives at each requested output
for ( plotnum = start_sequence; plotnum <= final_sequence;
plotnum = plotnum + step_sequence ) {
if ( deriv_x or deriv_y or deriv_z ) {
// loop over each field
for ( int var_num = 0; var_num <= fields.size()-1; var_num++ ) {
// check if field is a derivative field
bool input_deriv = false; // assume it's not a derivative
int var_len = fields[var_num].length();
if ( var_len > 2 ) {
if ( fields[var_num].substr(var_len-2,1) == "_" ) {
// // if second last char is an underscore then its a derivative field
input_deriv = true;
}
}
// parse for expansion type
find_expansion(grid_type, expan, fields[var_num]);
// read the field and setup for derivative
if ( fields[var_num] == "v" ) {
saved_v = true;
init_tracer_restart(fields[var_num],v);
gradient_op->setup_array(&v,expan[x_ind],expan[y_ind],expan[z_ind]); }
else if ( fields[var_num] == "w" ) {
saved_w = true;
init_tracer_restart(fields[var_num],w);
gradient_op->setup_array(&w,expan[x_ind],expan[y_ind],expan[z_ind]); }
else {
// else use u to hold the field
init_tracer_restart(fields[var_num],u);
gradient_op->setup_array(&u,expan[x_ind],expan[y_ind],expan[z_ind]);
}
if (master()) {
fprintf(stdout,"Expansions are (x,y,z): (%s, %s, %s)\n",
S_EXP_NAME[expan[x_ind]],S_EXP_NAME[expan[y_ind]],
S_EXP_NAME[expan[z_ind]]);
}
// X derivative
if (deriv_x) {
gradient_op->get_dx(&deriv_var,false);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max x derivative: %.6g\n",max_var);
// save the derivative
if ( !input_deriv ) {
snprintf(filename,100,"%s_x",fields[var_num].c_str()); }
else {
snprintf(filename,100,"%sx",fields[var_num].c_str());
}
write_array(deriv_var,filename,plotnum);
if (master())
fprintf(stdout,"Completed the write for %s.%d\n",filename,plotnum);
}
// Y derivative
if (deriv_y) {
gradient_op->get_dy(&deriv_var,false);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max y derivative: %.6g\n",max_var);
// save the derivative
if ( !input_deriv ) {
snprintf(filename,100,"%s_y",fields[var_num].c_str()); }
else {
snprintf(filename,100,"%sy",fields[var_num].c_str());
}
write_array(deriv_var,filename,plotnum);
if (master())
fprintf(stdout,"Completed the write for %s.%d\n",filename,plotnum);
}
// Z derivative
if (deriv_z) {
gradient_op->get_dz(&deriv_var,false);
double max_var = psmax(max(abs(deriv_var)));
if (master()) fprintf(stdout,"Max z derivative: %.6g\n",max_var);
// save the derivative
if ( !input_deriv ) {
snprintf(filename,100,"%s_z",fields[var_num].c_str()); }
else {
snprintf(filename,100,"%sz",fields[var_num].c_str());
}
write_array(deriv_var,filename,plotnum);
if (master())
fprintf(stdout,"Completed the write for %s.%d\n",filename,plotnum);
}
}
}
// 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) {
// u
init_tracer_restart("u",u);
// v
if ( !saved_v ) {
if ( v_exist ) {
init_tracer_restart("v",v); }
else {
if (master()) fprintf(stdout,"No v field, setting v=0\n");
v = 0;
}
}
// w
if ( !saved_w ) {
init_tracer_restart("w",w);
}
}
// X-component of vorticity
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);
write_array(deriv_var,"vortx",plotnum);
if (master())
fprintf(stdout,"Completed the write for vortx.%d\n",plotnum);
}
// Y-component of vorticity
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);
write_array(deriv_var,"vorty",plotnum);
if (master())
fprintf(stdout,"Completed the write for vorty.%d\n",plotnum);
}
// Z-component of vorticity
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);
write_array(deriv_var,"vortz",plotnum);
if (master())
fprintf(stdout,"Completed the write for vortz.%d\n",plotnum);
}
// Calculate Enstrophy
/*if ( do_enstrophy ) {
enstrophy_density(deriv_var, u, v, w, gradient_op, grid_type,
Nx, Ny, Nz);
double tot_enst = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*
(*get_quad_z())(kk)*deriv_var));
if (master()) fprintf(stdout,"Total Enstrophy: %.6g\n",tot_enst);
write_array(deriv_var,"enst",plotnum);
if (master())
fprintf(stdout,"Completed the write for enst.%d\n",plotnum);
}
// Calculate Viscous dissipation
if ( do_dissipation ) {
double mu = visco*rho_0; // dynamic viscosity
dissipation(deriv_var, u, v, w, gradient_op, grid_type,
Nx, Ny, Nz, mu);
double tot_diss = pssum(sum(