Implement surface stress calculation

The stress along the bottom and top surfaces are calculated
and written to a file (stresses_(top/bottom).txt) using the function
stresses_top and stresses_bottom in BaseCase. These functions are
especially useful with mapped cases.

src/BaseCase.cpp

 #include "BaseCase.hpp"
+#include "Science.hpp"
 #include "NSIntegrator.hpp"
 #include "TArray.hpp"
 #include <blitz/array.h>
@@ -452,3 +453,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==1 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==1 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);
+    }
+}

src/BaseCase.hpp

@@ -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;
@@ -170,6 +171,14 @@ class BaseCase {
       // 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

src/Science.cpp

@@ -426,7 +426,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];
         }
@@ -668,3 +668,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);
+    }
+}

src/Science.hpp

@@ -66,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)
 

src/cases/gravity_current/gravity_current.cpp

@@ -6,7 +6,6 @@
 // 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;
@@ -63,6 +62,8 @@ 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
 
@@ -191,10 +192,11 @@ class userControl : public BaseCase {
             iter++;
             // Set-up
             if ( iter == 1 ) {
-                if (compute_enstrophy or compute_dissipation) {
+                if ( compute_enstrophy or compute_dissipation or
+                        compute_stresses_top or compute_stresses_bottom ) {
                     temp1 = alloc_array(Nx,Ny,Nz);
                 }
-                // Determine last plot if restarting from the dump case
+                // Determine last plot if restarting from the dump file
                 if (restart_from_dump) {
                     next_plot = (restart_sequence+1)*plot_interval;    
                 }
@@ -332,6 +334,19 @@ class userControl : public BaseCase {
                         max_u,max_v,max_w,max_rho);
             }
 
+            // 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);
+            }
+            // 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);
+            }
+
             /* Write to disk if at correct time */
             if ((time - next_plot) > -1e-6) {
                 plot_number++;
@@ -446,7 +461,7 @@ int main(int argc, char ** argv) {
     option_category("Temporal options");
     add_option("final_time",&final_time,"Final time");
     add_option("plot_interval",&plot_interval,"Time between writes");
-    add_option("dt_max",&dt_max,0,"Maximum time step. Zero value results in the default");
+    add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
 
     option_category("Restart options");
     add_option("restart",&restarting,false,"Restart from prior output time.");
@@ -464,6 +479,8 @@ 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_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
     add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
 
     option_category("Filter options");
@@ -509,7 +526,7 @@ int main(int argc, char ** argv) {
     // Maximum buoyancy frequency (squared) if the initial stratification was stable
     N2_max = g*delta_rho/(2*delta_x);
     // Maximum time step
-    if (dt_max <= 0) {
+    if (dt_max == 0.0) {
         // if dt_max not given in spins.conf, use the buoyancy frequency
         dt_max = 0.5/sqrt(N2_max);
     }