Commit 6c338dbc authored by David Deepwell's avatar David Deepwell

Extend one diagnostic to handle mapped grids, and other minor fixes

The energy transfer rate from internal energy to BPE is generalized
to handle both mapped, and unmapped cases. What had been there
was in fact correct for both, though it wasn't realized at first.

Archived old wave_reader cases.

Minor fixes in other case files.
parent ae629429
......@@ -466,41 +466,23 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DT
// Internal energy converted to Background potential energy
// See Winters et al. 1995 (Available potential energy and mixing in density-stratified fluids)
// phi_i = - kappa * g * (surface integral of density at top boundary
// minus the surface integral of density at bottom boundary)
// phi_i = - kappa * g * (integral of density at top boundary
// minus the integral of density at bottom boundary)
void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
double kappa_rho, double rho_0, double g, int Nz,
bool dimensional_rho, bool mapped, TArrayn::DTArray * Hprime) {
double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho) {
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::Range all = blitz::Range::all();
if ( !mapped ) {
// Unmapped case
if ( !dimensional_rho ) {
phi_i = -kappa_rho * g * pssum(sum(
rho_0 * (rho(all,all,Nz-1) - rho(all,all,0))
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
} else {
phi_i = -kappa_rho * g * pssum(sum(
(rho(all,all,Nz-1) - rho(all,all,0))
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
}
if ( !dimensional_rho ) {
phi_i = -kappa_rho * g * pssum(sum(
rho_0 * (rho(all,all,Nz-1) - rho(all,all,0))
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
} else {
fprintf(stderr,"BPE from internal energy calculation is not available in mapped cases.\n");
assert(0);
// Mapped case
// Need to fix
/*if ( !dimensional_rho ) {
phi_i = -kappa_rho * g * pssum(sum(
rho_0 * ( rho(all,all,Nz-1) - rho(all,all,0)*pow(1+pow(*Hprime,2),0.5) )
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
} else {
phi_i = -kappa_rho * g * pssum(sum(
( rho(all,all,Nz-1) - rho(all,all,0)*pow(1+pow(*Hprime,2),0.5) )
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
}*/
phi_i = -kappa_rho * g * pssum(sum(
(rho(all,all,Nz-1) - rho(all,all,0))
* (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
}
}
......
......@@ -49,8 +49,7 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DT
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,
double kappa_rho, double rho_0, double g, int Nz,
bool dimensional_rho = false, bool mapped = false, TArrayn::DTArray * Hprime = NULL);
double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho = false);
// Quadrature weights
void compute_quadweights(int szx, int szy, int szz,
......
......@@ -124,7 +124,7 @@ class userControl : public BaseCase {
} else if (restarting and restart_from_dump) {
init_vels_dump(u, v, w);
} else{
// initial dipole (use v as the vector r^2)
// initial dipole
u = 0.5*U0 *
( (zz(kk)-Z1) * exp(-(pow(xx(ii)-X1,2) + pow(zz(kk)-Z1,2))/(r0*r0))
-(zz(kk)-Z2) * exp(-(pow(xx(ii)-X2,2) + pow(zz(kk)-Z2,2))/(r0*r0))
......@@ -295,7 +295,7 @@ class userControl : public BaseCase {
}
// Write to file
if (!(restarting and (iter==0)))
if (!(restarting and iter==0))
write_diagnostics(header, line, iter, restarting);
// and to the log file
fprintf(stdout,"[%d] (%.4g) %.4f: "
......
......@@ -154,7 +154,7 @@ class userControl : public BaseCase {
// Write the arrays and matlab readers
write_array(xg,"xgrid");
write_reader(xg,"xgrid",false);
if (Ny > 1 || rot_f != 0) {
if (Ny > 1) {
write_array(yg,"ygrid");
write_reader(yg,"ygrid",false);
}
......@@ -291,12 +291,7 @@ class userControl : public BaseCase {
// 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
}
compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
}
// viscous dissipation
double diss_tot = 0;
......@@ -385,7 +380,8 @@ class userControl : public BaseCase {
}
// Write to file
write_diagnostics(header, line, iter, restarting);
if (!(restarting and iter==0))
write_diagnostics(header, line, iter, restarting);
// and to the log file
fprintf(stdout,"[%d] (%.4g) %.4f: "
"%.4g %.4g %.4g %.4g\n",
......
......@@ -336,7 +336,7 @@ class userControl : public BaseCase {
}
// Write to file
if (!(restarting and (iter==0)))
if (!(restarting and iter==0))
write_diagnostics(header, line, iter, restarting);
// and to the log file
fprintf(stdout,"[%d] (%.4g) %.4f: "
......
......@@ -334,12 +334,7 @@ class userControl : public BaseCase {
// 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
}
compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
}
// viscous dissipation
double diss_tot = 0;
......@@ -428,7 +423,8 @@ class userControl : public BaseCase {
}
// Write to file
write_diagnostics(header, line, iter, restarting);
if (!(restarting and iter==0))
write_diagnostics(header, line, iter, restarting);
// and to the log file
fprintf(stdout,"[%d] (%.4g) %.4f: "
"%.4g %.4g %.4g %.4g\n",
......
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