Commit d21912d5 authored by Christopher Subich's avatar Christopher Subich
Browse files

Merge remote-tracking branch 'david/master'

parents 8ae27fe5 d317ebc8
#include "BaseCase.hpp"
#include "Science.hpp"
#include "NSIntegrator.hpp"
#include "TArray.hpp"
#include <blitz/array.h>
......@@ -412,7 +413,7 @@ void BaseCase::write_diagnostics(string header, string line,
FILE * diagnos_file = fopen("diagnostics.txt","a");
assert(diagnos_file);
// print header
if (iter == 1 and !restarting) {
if (iter == 0 and !restarting) {
fprintf(diagnos_file,"%s\n",clean_header.c_str());
}
// print the line of values
......@@ -421,6 +422,25 @@ void BaseCase::write_diagnostics(string header, string line,
fclose(diagnos_file);
}
// Write plot time information
void BaseCase::write_plot_times(double time, double clock_time, double comp_duration,
double avg_write_time, int plot_number, bool restarting) {
if (master()) {
// in log file
fprintf(stdout,"*Write time: %.6g. Average write time: %.6g.\n",
clock_time - comp_duration, avg_write_time);
// track in a file
FILE * plottimes_file = fopen("plot_times.txt","a");
assert(plottimes_file);
if ( plot_number==get_restart_sequence()+1 and !restarting )
fprintf(plottimes_file,"Output number, Simulation time (s), "
"Write time (s), Average write time (s)\n");
fprintf(plottimes_file,"%d, %.17f, %.12g, %.12g\n",
plot_number, time, clock_time - comp_duration, avg_write_time);
fclose(plottimes_file);
}
}
// 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) {
......@@ -452,3 +472,93 @@ void parse_boundary_conditions(const string xgrid_type, const string ygrid_type,
}
}
// Top stresses
void BaseCase::stresses_top(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const double mu, double time, int itercount, bool restarting) {
// set-up
static DTArray *tx = alloc_array(size_x(),size_y(),1);
static DTArray *ty = alloc_array(size_x(),size_y(),1);
blitz::firstIndex ii;
blitz::secondIndex jj;
// top stress ( along channel - x )
top_stress_x(*tx, u, temp, gradient_op, grid_type, size_z(), mu);
double top_tx_tot = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*(*tx)));
double top_tx_abs = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*abs(*tx)));
// top stress ( across channel - y )
top_stress_y(*ty, v, temp, gradient_op, grid_type, size_z(), mu);
double top_ty_tot = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*(*ty)));
double top_ty_abs = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*abs(*ty)));
// total top stress
double top_ts = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*pow(pow(*tx,2)+pow(*ty,2),0.5)));
// write to a stress diagnostic file
if (master()) {
FILE * stresses_file = fopen("stresses_top.txt","a");
assert(stresses_file);
if ( itercount==0 and !restarting )
fprintf(stresses_file,"Time, "
"Top_tx_tot, Top_tx_abs, Top_ty_tot, Top_ty_abs, Top_ts\n");
fprintf(stresses_file,"%.17f, "
"%.17g, %.17g, %.17g, %.17g, %.17g\n",
time,
top_tx_tot, top_tx_abs, top_ty_tot, top_ty_abs, top_ts);
fclose(stresses_file);
}
}
// Bottom stresses
void BaseCase::stresses_bottom(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const double mu, double time, int itercount, bool restarting) {
// set-up
static DTArray *tx = alloc_array(size_x(),size_y(),1);
static DTArray *ty = alloc_array(size_x(),size_y(),1);
blitz::firstIndex ii;
blitz::secondIndex jj;
// bottom stress ( along channel - x )
bottom_stress_x(*tx, Hprime, u, w, temp, gradient_op, grid_type, is_mapped(), mu);
double bot_tx_tot = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*(*tx)));
double bot_tx_abs = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*abs(*tx)));
// bottom stress ( across channel - y )
bottom_stress_y(*ty, Hprime, v, temp, gradient_op, grid_type, is_mapped(), mu);
double bot_ty_tot = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*(*ty)));
double bot_ty_abs = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*abs(*ty)));
// total bottom stress
double bot_ts = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*pow(pow(*tx,2)+pow(*ty,2),0.5)));
// write to a stress diagnostic file
if (master()) {
FILE * stresses_file = fopen("stresses_bottom.txt","a");
assert(stresses_file);
if ( itercount==0 and !restarting )
fprintf(stresses_file,"Time, "
"Bottom_tx_tot, Bottom_tx_abs, Bottom_ty_tot, Bottom_ty_abs, Bottom_ts\n");
fprintf(stresses_file,"%.17f, "
"%.17g, %.17g, %.17g, %.17g, %.17g\n",
time,
bot_tx_tot, bot_tx_abs, bot_ty_tot, bot_ty_abs, bot_ts);
fclose(stresses_file);
}
}
......@@ -6,6 +6,7 @@
#include <blitz/array.h>
#include "TArray.hpp"
#include "NSIntegrator.hpp"
#include "Science.hpp" // Science content
using namespace TArrayn;
using namespace NSIntegrator;
......@@ -166,10 +167,20 @@ class BaseCase {
template <class T> void add_diagnostic(const string str, const T val,
string & header, string & line);
void write_diagnostics(string header, string line, int iter, bool restarting);
void write_plot_times(double time, double clock_time, double comp_duration,
double avg_write_time, int plot_number, bool restarting);
// Generate an automatic grid for unmapped cases
virtual void automatic_grid(double MinX, double MinY, double MinZ,
Array<double,1> *xx=0, Array<double,1> *yy=0, Array<double,1> *zz=0);
// Surface Stresses
void stresses_top(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const double mu, double time, int itercount, bool restarting);
void stresses_bottom(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const double mu, double time, int itercount, bool restarting);
};
#include "BaseCase_impl.cc" // Include the implementation of the add_diagnostic template
......
......@@ -5,7 +5,7 @@
void WriteCaseFileSource(void)
{
char* filename;
if ( strcmp(casefilename, "cases/derivatives/derivatives.cpp") ) {
if ( strcmp(casefilename, "cases/derivatives/derivatives.cpp") == 0 ) {
filename = "derivatives.cpp";
} else {
filename = "spinscase.cpp";
......
......@@ -35,6 +35,9 @@ namespace NSIntegrator {
/* Control framework for a full run. At the first timestep(s), take
fractional timesteps to start up. */
// Analyze the initial conditions
usercode->analysis(times[0],us[0],vs[0],ws[0],tracers_now,pressure);
while (times[0] < (fintime - 1e-8*fabs(fintime))) {
// Since we're not done, we need to take a timestep.
......@@ -70,9 +73,6 @@ namespace NSIntegrator {
double desttime = times[0] + the_timestep;
// Analyze the initial conditions
usercode->analysis(times[0],us[0],vs[0],ws[0],tracers_now,pressure);
/* Now, take timesteps until we reach our one-timestep destination*/
bool starting_step = false;
......
......@@ -222,6 +222,32 @@ void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTA
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);
}
// Enstrophy Density: 1/2*(vort_x^2 + vort_y^2 + vort_z^2)
void enstrophy_density(TArrayn::DTArray & enst, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
const int Nx, const int Ny, const int Nz) {
// initalize temporary array
static DTArray *temp = alloc_array(Nx,Ny,Nz);
// square vorticity components
compute_vort_x(v, w, *temp, gradient_op, grid_type);
enst = pow(*temp,2);
compute_vort_y(u, w, *temp, gradient_op, grid_type);
enst += pow(*temp,2);
compute_vort_z(u, v, *temp, gradient_op, grid_type);
enst += pow(*temp,2);
enst = 0.5*enst;
}
// Viscous dissipation: 2*mu*e_ij*e_ij
void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
......@@ -280,21 +306,12 @@ void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray
diss *= 2.0*visco;
}
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);
}
bool compare_pairs( pair<double, double> a, pair<double, double> b ) {
return a.first < b.first;
}
// Compute Background Potential Energy (BPE)
void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DTArray & quad3,
int Nx, int Ny, int Nz, double Lx, double Ly, double Lz, double g,
double rho_0, int iter, bool dimensional_rho, bool mapped, Array<double,1> hill) {
// Tensor variables for indexing
......@@ -308,7 +325,6 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
// Arrays for sorting
static Array<double,3> sort_rho, sort_quad;
static DTArray *quad3; // quadrature weights
static vector<double> sort_hill(Nx), sort_dx(Nx);
static vector<double> Lx_partsum(Nx), hillvol_partsum(Nx);
double *local_vols = new double[numprocs];
......@@ -318,15 +334,9 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
static vector < pair<double, double> > height_width(Nx);
// Stuff to do once at the beginning
if (iter == 1) {
// create array of voxels
quad3 = alloc_array(Nx,Ny,Nz);
*quad3 = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
if (iter == 0) {
// adjust if mapped
if ( mapped ) {
*quad3 = (*quad3)*(Lz-hill(ii))/Lz;
// information about the hill
hill_vol = pssum(sum(hill*Ly*(*get_quad_x())));
hill_max = psmax(max(hill));
......@@ -379,7 +389,7 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
}
// Compute sorted rho
sortarray(rho, *quad3, sort_rho, sort_quad);
sortarray(rho, quad3, sort_rho, sort_quad);
// Volume of memory-local portion of sorted array
double local_vol = sum(sort_quad);
......@@ -409,7 +419,7 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
while ( (base_vol > Lx_partsum[II]*sort_hill[II] - hillvol_partsum[II]) and (II<Nx-1) ) {
II++;
}
assert(II>0 && "Something bad happened leading up to the tmpH calculation.\n");
assert(II>0 && "Something bad happened leading up to the tmpH calculation.");
// now subtract off the bit we went over by (draw yourself a picture, it works)
tmpH = sort_hill[II] - (sort_hill[II]*Lx_partsum[II] - base_vol - hillvol_partsum[II])/Lx_partsum[II-1];
}
......@@ -651,3 +661,149 @@ S_EXP swap_trig( S_EXP the_exp ) {
MPI_Finalize(); exit(1); // stop
}
}
// Bottom slope
void bottom_slope(TArrayn::DTArray & Hprime, TArrayn::DTArray & zgrid,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nx, const int Ny, const int Nz) {
// Set-up
DTArray & z_x = *alloc_array(Nx,Ny,Nz);
blitz::Range all = blitz::Range::all();
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::thirdIndex kk;
S_EXP expan[3];
assert(gradient_op);
// get bottom topography
Array<double,1> xx(split_range(Nx));
xx = zgrid(all,0,0);
// put into temp array, and take derivative
temp = xx(ii) + 0*jj + 0*kk;
find_expansion(grid_type, expan, "zgrid");
gradient_op->setup_array(&temp,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&z_x);
// flatten to get 2D array
Hprime(all,all,0) = z_x(all,all,0);
delete &z_x, xx;
}
// Top Stress (along channel - x)
void top_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & u,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nz, const double visco) {
// Set-up
blitz::Range all = blitz::Range::all();
S_EXP expan[3];
assert(gradient_op);
assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
// du/dz
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,false);
// top stress
stress_x(all,all,0) = -visco*temp(all,all,Nz-1);
}
// Top Stress (across channel - y)
void top_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & v,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nz, const double visco) {
// Set-up
blitz::Range all = blitz::Range::all();
S_EXP expan[3];
assert(gradient_op);
assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
// dv/dz
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,false);
// top stress
stress_y(all,all,0) = -visco*temp(all,all,Nz-1);
}
// Bottom Stress (along channel - x)
void bottom_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & Hprime,
TArrayn::DTArray & u, TArrayn::DTArray & w, TArrayn::DTArray & temp,
TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
const double visco) {
// Set-up
blitz::Range all = blitz::Range::all();
S_EXP expan[3];
assert(gradient_op);
assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
// Along channel bottom stress can also be slightly inaccurate when x boundary condition is free slip
// since u or w do not necessarily satisfy Neumann BCs and the derivative has Gibbs at domain edges
// This is only true if there is velocity at left or right (in x) end wall corner, so is likely always small
if (mapped) {
// -du/dx
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp,false);
temp = (-1)*temp;
// dw/dz
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,true);
// 2H'*(w_z-u_x)
stress_x(all,all,0) = 2*Hprime(all,all,0)*temp(all,all,0);
// dw/dx
find_expansion(grid_type, expan, "w");
gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp,false);
// du/dz
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,true);
// (1-(H')^2)*(u_z+w_x)
stress_x(all,all,0) += (1-pow(Hprime(all,all,0),2))*temp(all,all,0);
// multiply by mu/(1+(H')^2)
stress_x = visco/(1+pow(Hprime,2))*stress_x;
} else {
// du/dz
find_expansion(grid_type, expan, "u");
gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,false);
// bottom stress
stress_x(all,all,0) = visco*temp(all,all,0);
}
}
// Bottom Stress (across channel - y)
void bottom_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & Hprime,
TArrayn::DTArray & v, TArrayn::DTArray & temp,
TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
const double visco) {
// Set-up
blitz::Range all = blitz::Range::all();
S_EXP expan[3];
assert(gradient_op);
assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
// Across channel bottom stress can also be slightly inaccurate when x boundary condition is free slip
// since v does not necessarily satisfy Neumann BCs and the derivative has Gibbs at domain edges
// This is only true if there is velocity at left or right (in x) end wall corner, so is likely always small
if (mapped) {
// dv/dx
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dx(&temp,false);
// -v_x*H'
stress_y(all,all,0) = -temp(all,all,0)*Hprime(all,all,0);
// dv/dz
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,false);
// add to -v_x*H'
stress_y(all,all,0) = temp(all,all,0) + stress_y(all,all,0);
// multiply by mu/sqrt(1+(H')^2)
stress_y = visco/pow(1+pow(Hprime,2),0.5)*stress_y;
} else {
// dv/dz
find_expansion(grid_type, expan, "v");
gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
gradient_op->get_dz(&temp,false);
// bottom stress
stress_y(all,all,0) = visco*temp(all,all,0);
}
}
......@@ -34,14 +34,18 @@ void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTA
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);
// 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,
const int Nx, const int Ny, const int Nz);
// Viscous dissipation
void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
const int Nx, const int Ny, const int Nz, const double visco);
// Background Potential Energy (BPE)
void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, int Nx, int Ny, int Nz,
double Lx, double Ly, double Lz, double g, double rho_0, int iter,
void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DTArray & quad3,
int Nx, int Ny, int Nz, double Lx, double Ly, double Lz, double g, double rho_0, int iter,
bool dimensional_rho = false, bool mapped = false, Array<double,1> hill = Array<double,1>());
// Internal energy converted to BPE
void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
......@@ -62,6 +66,27 @@ void find_expansion(const string * grid_type, Transformer::S_EXP * expan, string
// switch trig function
Transformer::S_EXP swap_trig( Transformer::S_EXP the_exp );
// Bottom slope
void bottom_slope(TArrayn::DTArray & Hprime, TArrayn::DTArray & zgrid,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nx, const int Ny, const int Nz);
// Top stresses
void top_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & u,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nz, const double visco);
void top_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & v,
TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
const string * grid_type, const int Nz, const double visco);
// Bottom stresses
void bottom_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & Hprime,
TArrayn::DTArray & u, TArrayn::DTArray & w, TArrayn::DTArray & temp,
TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
const double visco);
void bottom_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & Hprime,
TArrayn::DTArray & v, TArrayn::DTArray & temp,
TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
const double visco);
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
......
......@@ -243,7 +243,7 @@ class userControl : public BaseCase {
}
// Calculate Enstrophy
/*if ( do_enstrophy ) {
if ( do_enstrophy ) {
enstrophy_density(deriv_var, u, v, w, gradient_op, grid_type,
Nx, Ny, Nz);
double tot_enst = pssum(sum(
......@@ -255,6 +255,7 @@ class userControl : public BaseCase {
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
......@@ -268,7 +269,7 @@ class userControl : public BaseCase {
write_array(deriv_var,"diss",plotnum);
if (master())
fprintf(stdout,"Completed the write for diss.%d\n",plotnum);
}*/
}
}
}
......
/* Script for the formation of a gravity current with zero initial velocity
* and no topography */
/* Script for the formation of a gravity current with:
* zero initial velocity
* density stratification (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 "../../Science.hpp" // Science content
#include <random/normal.h> // Blitz random number generator
using namespace ranlib;
......@@ -22,6 +23,8 @@ blitz::thirdIndex kk;
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];
......@@ -40,6 +43,11 @@ double delta_rho; // density difference between different layers (
double delta_x; // horizontal transition length (m)
double Lmix; // Width of mixed region (m)
// Topography parameters
double hill_height; // height of hill (m)
double hill_centre; // position of hill peak (m)
double hill_width; // width of hill (m)
// Temporal parameters
double final_time; // Final time (s)
double plot_interval; // Time between field writes (s)
......@@ -63,6 +71,9 @@ 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
// Maximum squared buoyancy frequency
......@@ -72,12 +83,15 @@ double N2_max;
class userControl : public BaseCase {
public:
// Grid arrays
Array<double,1> xx, yy, zz;
// Grid and topography arrays
DTArray *xgrid, *ygrid, *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
// Arrays and operators for derivatives
Grad * gradient_op;
DTArray *temp1;
DTArray *temp1, *dxdydz;
// Timing variables (for outputs and measuring time steps)
int plot_number; // plot output number
......@@ -119,6 +133,39 @@ class userControl : public BaseCase {
// Number of tracers (the first is density)
int numtracers() const { return Num_tracers; }
// Create mapped grid
bool is_mapped() const { return mapped; }
void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
xgrid = alloc_array(Nx,Ny,Nz);
ygrid = alloc_array(Nx,Ny,Nz);
zgrid = alloc_array(Nx,Ny,Nz);
// over-write zz to be between -1 and 1
// (zz defined in automatic_grid below)
zz = -cos(ii*M_PI/(Nz-1)); // Chebyshev in vertical