Commit d317ebc8 authored by David Deepwell's avatar David Deepwell
Browse files

Update gravity_current.cpp to allow mapping

parent 17c63523
......@@ -311,7 +311,7 @@ bool compare_pairs( pair<double, double> a, pair<double, double> b ) {
}
// 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
......@@ -325,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];
......@@ -336,14 +335,8 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
// Stuff to do once at the beginning
if (iter == 0) {
// create array of voxels
quad3 = alloc_array(Nx,Ny,Nz);
*quad3 = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
// 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));
......@@ -396,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);
......
......@@ -44,8 +44,8 @@ void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray
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,
......
/* 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 --------------------- */
......@@ -21,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];
......@@ -39,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)
......@@ -74,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
Grad * gradient_op;
DTArray *temp1;
DTArray *temp1, *dxdydz;
// Timing variables (for outputs and measuring time steps)
int plot_number; // plot output number
......@@ -121,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
// Define topography
topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
// make full grids
xg = xx(ii) + 0*jj + 0*kk;
yg = yy(jj) + 0*ii + 0*kk;
zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
*xgrid = xg;
*ygrid = yg;
*zgrid = zg;
// Write the arrays and matlab readers
write_array(xg,"xgrid");
write_reader(xg,"xgrid",false);
if (Ny > 1 || rot_f != 0) {
write_array(yg,"ygrid");
write_reader(yg,"ygrid",false);
}
write_array(zg,"zgrid");
write_reader(zg,"zgrid",false);
}
/* Initialize velocities */
void init_vels(DTArray & u, DTArray & v, DTArray & w) {
if (master()) fprintf(stdout,"Initializing velocities\n");
......@@ -170,6 +215,8 @@ class userControl : public BaseCase {
} else {
// Density configuration
rho = delta_rho*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
// Important! if mapped, and rho depends on z
// then (*zgrid)(ii,jj,kk), must be used in stead of zz(kk)
// Write the array
write_array(rho,"rho",plot_number);
}
......@@ -194,9 +241,25 @@ class userControl : public BaseCase {
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;
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
......@@ -210,27 +273,34 @@ class userControl : public BaseCase {
// 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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
}
if ( Nz > 1 ) {
ke_z = pssum(sum(0.5*rho_0*(w*w)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
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 pe_tot = pssum(sum(rho_0*(1+*tracers[RHO])*g*(zz(kk) - MinZ)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
double BPE_tot = 0;
if (compute_BPE) {
compute_Background_PE(BPE_tot, *tracers[RHO], Nx, Ny, Nz, Lx, Ly, Lz, g, rho_0, iter);
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) {
compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
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;
......@@ -238,8 +308,7 @@ class userControl : public BaseCase {
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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
diss_tot = pssum(sum((*temp1)*(*dxdydz)));
}
// Vorticity / Enstrophy
double max_vort_x = 0, enst_x_tot = 0;
......@@ -250,22 +319,19 @@ class userControl : public BaseCase {
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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
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)*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
}
}
// max of fields
......@@ -275,8 +341,7 @@ class userControl : public BaseCase {
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]))*
(*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
double mass = pssum(sum(rho_0*(1+(*tracers[RHO]))*(*dxdydz)));
if (master()) {
// add diagnostics to buffers
......@@ -334,22 +399,18 @@ class userControl : public BaseCase {
// Top Surface Stresses
if ( compute_stresses_top ) {
static DTArray & Hprime = *alloc_array(Nx,Ny,1);
Hprime = 0*ii + 0*jj;
stresses_top(u, v, w, Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
}
// Bottom Surface Stresses
if ( compute_stresses_bottom ) {
static DTArray & Hprime = *alloc_array(Nx,Ny,1);
Hprime = 0*ii + 0*jj;
stresses_bottom(u, v, w, Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
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)
//Write the arrays
// Write the arrays
write_array(u,"u",plot_number);
write_array(w,"w",plot_number);
if (Ny > 1 || rot_f != 0)
......@@ -389,7 +450,7 @@ class userControl : public BaseCase {
// Constructor: Initialize local variables
userControl():
xx(split_range(Nx)), yy(Ny), zz(Nz),
gradient_op(0),
topo(split_range(Nx)), gradient_op(0),
plot_number(restart_sequence),
next_plot(initial_time + plot_interval)
{ compute_quadweights(
......@@ -432,6 +493,12 @@ int main(int argc, char ** argv) {
" NO_SLIP: Chebyhsev expansion");
add_option("type_y",&ygrid_type,"FOURIER","Grid type in Y");
add_option("type_z",&zgrid_type,"Grid type in Z");
add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
option_category("Topography parameters");
add_option("hill_height",&hill_height,"Height of hill");
add_option("hill_centre",&hill_centre,"location of hill peak");
add_option("hill_width",&hill_width,"Width of hill");
option_category("Physical parameters");
add_option("g",&g,9.81,"Gravitational acceleration");
......@@ -466,7 +533,7 @@ int main(int argc, char ** argv) {
add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
"Calculate BPE gained from internal energy?");
add_option("compute_stresses_top", &compute_stresses_top, false,"Calculate top surfaces stresses?");
add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
......
......@@ -3,7 +3,7 @@
# Spatial Parameters
Lx = 1.0
Ly = 0.1
Lz = 0.1
Lz = 0.5
Nx = 256
Ny = 1
Nz = 128
......@@ -29,9 +29,22 @@ delta_rho = 0.005
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 50
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 = 1e-3
......@@ -41,9 +54,10 @@ f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
#compute_BPE = true
#compute_internal_to_BPE = true
#compute_stresses_top = false
#compute_stresses_bottom = false
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