Create BPE calculator

The background Potential Energy (BPE) is calculated in
compute_Background_PE in Science.cpp. The APE can
be calculated afterwards from the BPE and PE from the
diagnostic file.

src/Science.cpp

@@ -9,6 +9,7 @@
 #include "Split_reader.hpp"
 #include "T_util.hpp"
 #include "Parformer.hpp"
+#include "Sorter.hpp"
 
 // Marek's Overturning Diagnostic
 
@@ -229,6 +230,67 @@ void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArra
     compute_vort_z(vortz, u, v, gradient_op, grid_type);
 }
 
+// Compute 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 g, double rho_0, int iter) {
+    // Tensor variables for indexing
+    blitz::firstIndex ii;
+    blitz::secondIndex jj;
+    blitz::thirdIndex kk;
+    // Set parameters for sorting
+    int myrank, numprocs;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+    double *local_vols;
+    double *local_vols_r;
+    local_vols   = new double[numprocs];
+    local_vols_r = new double[numprocs];
+
+    // Arrays for sorting
+    static Array<double,3> sort_rho, sort_quad;
+    static DTArray *quad3;     // quadrature weights
+    if (iter == 1) {
+        quad3 = alloc_array(Nx,Ny,Nz);
+        *quad3 = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+    }
+
+    // Compute sorted rho
+    sortarray(rho, *quad3, sort_rho, sort_quad);
+
+    // Volume of memory-local portion of sorted array
+    double local_vol = sum(sort_quad);
+
+    // Share local volume information with other processors
+    for (int II = 0; II < numprocs; II++) {
+        local_vols[II] = (II <= myrank) ? local_vol : 0;
+    }
+    MPI_Allreduce(local_vols, local_vols_r, numprocs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    // The sum of the volumes with higher density is local_vols_r[rank]
+    double base_vol = local_vols_r[myrank];
+
+    // Now compute BPE from the sorted rho field
+    double tmpH = base_vol;
+    double dH = 0;
+    double BPE = 0;
+    for (int II = sort_rho.lbound(firstDim); II <= sort_rho.ubound(firstDim); II++) {
+        for (int KK = sort_rho.lbound(thirdDim); KK <= sort_rho.ubound(thirdDim); KK++) {
+            for (int JJ = sort_rho.lbound(secondDim); JJ <= sort_rho.ubound(secondDim); JJ++) {
+                dH = sort_quad(II,JJ,KK)/(Lx*Ly);
+                BPE += g*rho_0*(1+sort_rho(II,JJ,KK))*(tmpH - 0.5*dH)*dH*Lx*Ly;
+                tmpH -= dH;
+            }
+        }
+    }
+
+    // Share BPE with other processors
+    MPI_Allreduce(&BPE, &BPE_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    // delete variables
+    delete[] local_vols;
+    delete[] local_vols_r;
+}
+
 // Global arrays to store quadrature weights
 Array<double,1> _quadw_x, _quadw_y, _quadw_z;
 

src/Science.hpp

@@ -35,6 +35,9 @@ void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArra
         TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
         TArrayn::Grad * gradient_op, const string * grid_type);
 
+// Background Potential Energy
+void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, int Nx, int Ny, int Nz,
+        double Lx, double Ly, double g, double rho_0, int iter);
 
 // Quadrature weights
 void compute_quadweights(int szx, int szy, int szz, 

src/cases/gravity_current/gravity_current.cpp

@@ -7,7 +7,7 @@
 #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
+#include <random/normal.h>         // Blitz random number generator
 
 using namespace ranlib;
 
@@ -60,6 +60,7 @@ double compute_start_time;      // real (clock) time when computation begins (af
 // other options
 double perturb;                 // Initial velocity perturbation
 bool compute_enstrophy;         // Compute enstrophy?
+bool compute_BPE;               // Compute enstrophy?
 int iter = 0;                   // Iteration counter
 
 // Maximum squared buoyancy frequency
@@ -211,6 +212,10 @@ class userControl : public BaseCase {
                         (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
             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;
+            if (compute_BPE) {
+                compute_Background_PE(BPE_tot, *tracers[RHO], Nx, Ny, Nz, Lx, Ly, g, rho_0, iter);
+            }
             // Vorticity / Enstrophy
             double max_vort_x, enst_x_tot;
             double max_vort_y, enst_y_tot;
@@ -253,6 +258,9 @@ class userControl : public BaseCase {
                 add_diagnostic("KE_x", ke_x,            header, line);
                 add_diagnostic("KE_z", ke_z,            header, line);
                 add_diagnostic("PE_tot", pe_tot,        header, line);
+                if (compute_BPE) {
+                    add_diagnostic("BPE_tot", BPE_tot,      header, line);
+                }
                 if (compute_enstrophy) {
                     add_diagnostic("Max_vort_y", max_vort_y,    header, line);
                     add_diagnostic("Enst_y_tot", enst_y_tot,    header, line);
@@ -403,6 +411,7 @@ int main(int argc, char ** argv) {
     option_category("Other options");
     add_option("perturb",&perturb,"Initial perturbation in velocity");
     add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+    add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
 
     option_category("Filter options");
     add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");