Commit 8e790761 authored by David Deepwell's avatar David Deepwell
Browse files

Update wave_reader case

parent d21912d5
......@@ -319,6 +319,27 @@ void BaseCase::init_tracer_dump(const std::string & field, DTArray & the_tracer)
return;
}
/* Read field from input data type */
void BaseCase::init_field(const std::string & field,
const std::string & filename, DTArray & the_field, input_types input_data_type) {
if (input_data_type == MATLAB) {
// from matlab
if (master()) fprintf(stdout,"Reading MATLAB-format %s (%d x %d) from %s\n",
field.c_str(),size_x(),size_z(),filename.c_str());
read_2d_slice(the_field,filename.c_str(),size_x(),size_z());
} else if (input_data_type == CTYPE) {
// from 2D spins output (generally, to go to 3D)
if (master()) fprintf(stdout,"Reading CTYPE-format %s (%d x %d) from %s\n",
field.c_str(),size_x(),size_z(),filename.c_str());
read_2d_restart(the_field,filename.c_str(),size_x(),size_z());
} else {
// from full (2D or 3D) spins output
// only suitable for initializing the grid
if (master()) fprintf(stdout,"Reading %s from %s\n",field.c_str(),filename.c_str());
read_array(the_field,filename.c_str(),size_x(),size_y(),size_z());
}
}
/* write out vertical chain of data */
void BaseCase::write_chain(const char *filename, DTArray & val, int Iout, int Jout, double time) {
FILE *fid=fopen(filename,"a");
......@@ -441,6 +462,18 @@ void BaseCase::write_plot_times(double time, double clock_time, double comp_dura
}
}
// parse file types
void parse_datatype(const string datatype, input_types & input_data_type) {
if (datatype=="MATLAB") { input_data_type = MATLAB; }
else if (datatype == "CTYPE") { input_data_type = CTYPE; }
else if (datatype == "FULL") { input_data_type = FULL; }
else {
if (master())
fprintf(stderr,"Invalid option %s received for file_type\n",datatype.c_str());
MPI_Finalize(); exit(1);
}
}
// parse expansion types
void parse_boundary_conditions(const string xgrid_type, const string ygrid_type,
const string zgrid_type, DIMTYPE & intype_x, DIMTYPE & intype_y, DIMTYPE & intype_z) {
......
......@@ -13,6 +13,13 @@ using namespace NSIntegrator;
using blitz::Array;
using std::vector;
// Possible input data types
static enum input_types {
MATLAB,
CTYPE,
FULL
} input_data_types;
class BaseCase {
/* To reduce boilerplate, wrap some of the long functions, only calling
them if actually used by usercode. For example, a tracer-free code
......@@ -98,6 +105,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_field(const std::string & field, const std::string & filename,
DTArray & the_field, input_types input_data_type);
virtual void init_grid_restart(const std::string & component,
const std::string & filename, DTArray & grid);
......@@ -191,6 +200,8 @@ extern template void BaseCase::add_diagnostic<int>(const string str, const int v
extern template void BaseCase::add_diagnostic<double>(const string str, const double val,
string & header, string & line);
// parse for data type
void parse_datatype(const string datatype, input_types & input_data_type);
// parse expansion types
void parse_boundary_conditions(const string xgrid_type, const string ygrid_type,
const string zgrid_type, DIMTYPE & intype_x, DIMTYPE & intype_y, DIMTYPE & intype_z);
......
## density wave reader configuration file
# Spatial Parameters
Lx = 1.0
Ly = 0.1
Lz = 0.5
Nx = 256
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_rho = 1e-6
# 2D Input Files
file_type = CTYPE
xgrid = x2d
zgrid = z2d
u_file = u2d
w_file = w2d
rho_file = rho2d
tracer_file = trc2d
# Tracer Parameters
enable_tracer = false
tracer_kappa = 1.0e-6
tracer_gravity = 0
# Temporal Parameters
final_time = 10
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 0e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
#compute_BPE = true
#compute_internal_to_BPE = true
#compute_stresses_top = false
#compute_stresses_bottom = false
/* wave_reader.cpp -- general case for
* extending 2D into 3D
* starting from Matlab initialized data
* or restarting any case
* density stratifications (not salt/temp)
* and with or without topography */
/* ------------------ Top matter --------------------- */
// Required headers
#include "../../BaseCase.hpp" // Support file containing default implementations of several functions
#include "../../Options.hpp" // config-file parser
#include <random/normal.h> // Blitz random number generator
using namespace ranlib;
// 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
double MinX, MinY, MinZ; // Minimum x/y/z points (m)
// Mapped grid?
bool mapped;
// Grid types
DIMTYPE intype_x, intype_y, intype_z;
string grid_type[3];
// Physical parameters
double g, rot_f, rho_0; // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
double visco; // viscosity (m^2/s)
double mu; // dynamic viscosity (kg/(m·s))
double kappa_rho; // diffusivity of density (m^2/s)
// tracer options
const int RHO = 0; // index for rho
const int TRCR = 1; // index for tracer
int Num_tracers; // number of tracers (density and dyes)
bool is_tracer; // is there a tracer?
double kappa_trc, tracer_g; // tracer diffusivity (m^2/s) and tracer gravitational acceleration
// Temporal parameters
double final_time; // Final time (s)
double plot_interval; // Time between field writes (s)
double dt_max; // maximum time step (s)
// Restarting options
bool restarting; // are you restarting?
double initial_time; // initial start time of simulation
int restart_sequence; // output number to restart from
// Dump parameters
bool restart_from_dump; // restarting from dump?
double compute_time; // requested computation time
double avg_write_time; // average time to write all output fields at one output
double real_start_time; // real (clock) time when simulation begins
double compute_start_time; // real (clock) time when computation begins (after initialization)
// other options
double perturb; // Initial velocity perturbation
bool compute_enstrophy; // Compute enstrophy?
bool compute_dissipation; // Compute dissipation?
bool compute_BPE; // Compute background potential energy?
bool compute_internal_to_BPE; // Compute BPE gained from internal energy?
bool compute_stresses_top; // Compute top surface stresses?
bool compute_stresses_bottom; // Compute bottom surface stresses?
bool write_pressure; // Write the pressure field?
int iter = 0; // Iteration counter
// Input file names
string xgrid_filename,
ygrid_filename,
zgrid_filename,
u_filename,
v_filename,
w_filename,
rho_filename,
tracer_filename;
/* ------------------ Adjust the class --------------------- */
class userControl : public BaseCase {
public:
// Grid and topography arrays
DTArray *zgrid; // Full grid fields
Array<double,1> xx, yy, zz; // 1D grid vectors
Array<double,1> topo; // topography vector
DTArray *Hprime; // derivative of topography vector
// Arrays and operators for derivatives
Grad * gradient_op;
DTArray *temp1, *dxdydz;
// Timing variables (for outputs and measuring time steps)
int plot_number; // plot output number
double next_plot; // time of next output write
double comp_duration; // clock time since computation began
double clock_time; // current clock time
// 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; }
// Coriolis parameter, viscosity, and diffusivities
double get_rot_f() const { return rot_f; }
double get_visco() const { return visco; }
double get_diffusivity(int t_num) const {
switch (t_num) {
case RHO:
return kappa_rho;
case TRCR:
return kappa_trc;
default:
if (master()) fprintf(stderr,"Invalid tracer number!\n");
MPI_Finalize(); exit(1);
}
}
// Temporal parameters
double init_time() const { return initial_time; }
int get_restart_sequence() const { return restart_sequence; }
double get_dt_max() const { return dt_max; }
double get_next_plot() { return next_plot; }
// Number of tracers (the first is density)
int numActive() const { return 1; }
int numPassive() const {
if (is_tracer) return 1;
else return 0;
}
// read grid if mapped
bool is_mapped() const { return mapped; }
void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
zgrid = alloc_array(Nx,Ny,Nz);
blitz::Range all = blitz::Range::all();
// read grid from file
init_field("xgrid", xgrid_filename, xg, input_data_types);
init_field("zgrid", zgrid_filename, zg, input_data_types);
*zgrid = zg;
// Automatically generate y-grid
yg = 0*ii + MinY + Ly*(0.5+jj)/Ny + 0*kk;
// Write the arrays and matlab readers
write_array(xg,"xgrid");
write_array(zg,"zgrid");
write_reader(xg,"xgrid",false);
write_reader(zg,"zgrid",false);
if ( Ny > 1 ) {
write_array(yg,"xgrid");
write_reader(yg,"xgrid",false);
}
// Define topography
topo = zg(all,0,0);
}
/* Initialize velocities */
void init_vels(DTArray & u, DTArray & v, DTArray & w) {
if (master()) fprintf(stdout,"Initializing velocities\n");
// if restarting
if (restarting and !restart_from_dump) {
init_vels_restart(u, v, w);
} else if (restarting and restart_from_dump) {
init_vels_dump(u, v, w);
} else {
// error check
if (input_data_types == FULL) {
if (master())
fprintf(stderr,"FULL option chosen, turn restart on\n");
MPI_Finalize(); exit(1);
}
// else start from other data formats
init_field("u", u_filename, u, input_data_types);
init_field("w", w_filename, w, input_data_types);
v = 0*ii + 0*jj + 0*kk;
// Add a random perturbation to trigger any 3D instabilities
int myrank;
MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
Normal<double> rnd(0,1);
for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
rnd.seed(i);
for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
u(i,j,k) *= 1+perturb*rnd.random();
w(i,j,k) *= 1+perturb*rnd.random();
if ( Ny > 1 || rot_f != 0)
v(i,j,k) *= 1+perturb*rnd.random();
}
}
}
// Write the arrays
write_array(u,"u",plot_number);
write_array(w,"w",plot_number);
if (Ny > 1 || rot_f != 0) {
write_array(v,"v",plot_number);
}
}
}
/* Initialize the tracers (density and dyes) */
void init_tracers(vector<DTArray *> & tracers) {
if (master()) fprintf(stdout,"Initializing tracers\n");
// Sanity checks
assert(numtracers() == int(tracers.size()));
assert(numtracers() >= 1);
// if restarting
if (restarting and !restart_from_dump) {
init_tracer_restart("rho",*tracers[RHO]);
if (is_tracer)
init_tracer_restart("tracer",*tracers[TRCR]);
} else if (restarting and restart_from_dump) {
init_tracer_dump("rho",*tracers[RHO]);
if (is_tracer)
init_tracer_dump("tracer",*tracers[TRCR]);
} else {
// error check
if (input_data_types == FULL) {
if (master())
fprintf(stderr,"FULL option chosen, turn restart on\n");
MPI_Finalize(); exit(1);
}
// else start from other data formats
init_field("rho", rho_filename, *tracers[RHO], input_data_types);
if (is_tracer)
init_field("tracer", tracer_filename, *tracers[TRCR], input_data_types);
// Write the arrays
write_array(*tracers[RHO],"rho",plot_number);
if (is_tracer)
write_array(*tracers[TRCR],"tracer",plot_number);
}
}
/* Forcing in the momentum equations */
void forcing(double t, const DTArray & u, DTArray & u_f,
const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
u_f = +rot_f*v;
v_f = -rot_f*u;
w_f = -g*(*tracers[RHO]); // tracers[RHO] is defined as rho/rho_0
*tracers_f[RHO] = 0;
if (is_tracer) {
*tracers_f[TRCR] = 0;
w_f += -tracer_g*(*tracers[TRCR]);
}
}
/* Basic analysis: compute secondary variables, and save fields and diagnostics */
void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
vector<DTArray *> & tracers, DTArray & pressure) {
// Set-up
if ( iter == 0 ) {
if ( compute_enstrophy or compute_dissipation or
compute_stresses_top or compute_stresses_bottom ) {
temp1 = alloc_array(Nx,Ny,Nz);
}
if ( compute_stresses_top or compute_stresses_bottom ) {
// initialize the vector of the bottom slope (Hprime)
Hprime = alloc_array(Nx,Ny,1);
if (is_mapped()) {
bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
} else {
topo = 0*ii;
*Hprime = 0*ii + 0*jj;
}
}
// Determine last plot if restarting from the dump file
if (restart_from_dump) {
next_plot = (restart_sequence+1)*plot_interval;
}
// initialize the size of each voxel
dxdydz = alloc_array(Nx,Ny,Nz);
*dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
if (is_mapped()) {
*dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
}
}
// update clocks
if (master()) {
clock_time = MPI_Wtime();
comp_duration = clock_time - compute_start_time;
}
/* Calculate and write out useful information */
// Energy (PE assumes density is density anomaly)
double ke_x = 0, ke_y = 0, ke_z = 0;
if ( Nx > 1 ) {
ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
}
if (Ny > 1 || rot_f != 0) {
ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
}
if ( Nz > 1 ) {
ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
}
double pe_tot;
if (is_mapped()) {
pe_tot = pssum(sum(rho_0*(1+(*tracers[RHO]))*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
} else {
pe_tot = pssum(sum(rho_0*(1+(*tracers[RHO]))*g*(zz(kk) - MinZ)*(*dxdydz)));
}
double BPE_tot = 0;
if (compute_BPE) {
compute_Background_PE(BPE_tot, *tracers[RHO], *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
g, rho_0, iter, false, is_mapped(), topo);
}
// Conversion from internal energy to background potential energy
double phi_i = 0;
if (compute_internal_to_BPE) {
if (!is_mapped()) {
compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz,
false, is_mapped(), Hprime);
} else {
// this is not finished yet for the mapped case
}
}
// viscous dissipation
double diss_tot = 0;
double max_diss = 0;
if (compute_dissipation) {
dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
max_diss = psmax(max(*temp1));
diss_tot = pssum(sum((*temp1)*(*dxdydz)));
}
// Vorticity / Enstrophy
double max_vort_x = 0, enst_x_tot = 0;
double max_vort_y = 0, enst_y_tot = 0;
double max_vort_z = 0, enst_z_tot = 0;
if (compute_enstrophy) {
// x-vorticity
if (Ny > 1 and Nz > 1) {
compute_vort_x(*temp1, v, w, gradient_op, grid_type);
max_vort_x = psmax(max(abs(*temp1)));
enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
}
// y-vorticity
if (Nx > 1 and Nz > 1) {
compute_vort_y(*temp1, u, w, gradient_op, grid_type);
max_vort_y = psmax(max(abs(*temp1)));
enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
}
// z-vorticity
if (Nx > 1 and Ny > 1) {
compute_vort_z(*temp1, u, v, gradient_op, grid_type);
max_vort_z = psmax(max(abs(*temp1)));
enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
}
}
// max of fields
double max_u = psmax(max(abs(u)));
double max_v = psmax(max(abs(v)));
double max_w = psmax(max(abs(w)));
double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
double max_rho = psmax(max(abs(*tracers[RHO])));
// total mass (tracers[RHO] is non-dimensional density)
double mass = pssum(sum(rho_0*(1+(*tracers[RHO]))*(*dxdydz)));
if (master()) {
// add diagnostics to buffers
string header, line;
add_diagnostic("Iter", iter, header, line);
add_diagnostic("Clock_time", comp_duration, header, line);
add_diagnostic("Time", time, header, line);
add_diagnostic("Max_vel", max_vel, header, line);
add_diagnostic("Max_density", max_rho, header, line);
add_diagnostic("Mass", mass, header, line);
add_diagnostic("PE_tot", pe_tot, header, line);
if (compute_BPE) {
add_diagnostic("BPE_tot", BPE_tot, header, line);
}
if (compute_internal_to_BPE) {
add_diagnostic("BPE_from_int", phi_i, header, line);
}
if (compute_dissipation) {
add_diagnostic("Max_diss", max_diss, header, line);
add_diagnostic("Diss_tot", diss_tot, header, line);
}
if (Nx > 1) {
add_diagnostic("Max_u", max_u, header, line);
add_diagnostic("KE_x", ke_x, header, line);
}
if (Ny > 1 || rot_f != 0) {
add_diagnostic("Max_v", max_v, header, line);
add_diagnostic("KE_y", ke_y, header, line);
}
if (Nz > 1) {
add_diagnostic("Max_w", max_w, header, line);
add_diagnostic("KE_z", ke_z, header, line);
}
if (Ny > 1 && Nz > 1 && compute_enstrophy) {
add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
add_diagnostic("Max_vort_x", max_vort_x, header, line);
}
if (Nx > 1 && Nz > 1 && compute_enstrophy) {
add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
add_diagnostic("Max_vort_y", max_vort_y, header, line);
}
if (Nx > 1 && Ny > 1 && compute_enstrophy) {
add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
add_diagnostic("Max_vort_z", max_vort_z, header, line);
}
// Write to file
write_diagnostics(header, line, iter, restarting);
// and to the log file
fprintf(stdout,"[%d] (%.4g) %.4f: "
"%.4g %.4g %.4g %.4g\n",
iter,comp_duration,time,
max_u,max_v,max_w,max_rho);
}
// Top Surface Stresses
if ( compute_stresses_top ) {
stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
}
// Bottom Surface Stresses
if ( compute_stresses_bottom ) {
stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
}
/* Write to disk if at correct time */
if ((time - next_plot) > -1e-6) {
plot_number++;
comp_duration = MPI_Wtime(); // time just before write (for dump)