Add viscous dissipation calculation

src/Science.cpp

@@ -222,6 +222,64 @@ void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTA
     gradient_op->get_dx(&vortz,true);
 }
 
+// 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,
+        const int Nx, const int Ny, const int Nz, const double visco) {
+    // Set-up
+    static DTArray *temp = alloc_array(Nx,Ny,Nz);
+    S_EXP expan[3];
+    assert(gradient_op);
+
+    // 1st term: e_11^2 = (du/dx)^2
+    find_expansion(grid_type, expan, "u");
+    gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+    gradient_op->get_dx(temp,false);
+    diss = pow(*temp,2);
+    // 2nd term: e_22^2 = (dv/dy)^2
+    find_expansion(grid_type, expan, "v");
+    gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+    gradient_op->get_dy(temp,false);
+    diss += pow(*temp,2);
+    // 3rd term: e_33^2 = (dw/dz)^2
+    find_expansion(grid_type, expan, "w");
+    gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+    gradient_op->get_dz(temp,false);
+    diss += pow(*temp,2);
+    // 4th term: 2e_12^2 = 2*(1/2*(u_y + v_x))^2
+    // u_y
+    find_expansion(grid_type, expan, "u");
+    gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+    gradient_op->get_dy(temp,false);
+    // v_x
+    find_expansion(grid_type, expan, "v");
+    gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+    gradient_op->get_dx(temp,true);
+    diss += 2.0*pow(0.5*(*temp),2);
+    // 5th term: 2e_13^2 = 2*(1/2*(u_z + w_x))^2
+    // u_z
+    find_expansion(grid_type, expan, "u");
+    gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+    gradient_op->get_dz(temp,false);
+    // w_x
+    find_expansion(grid_type, expan, "w");
+    gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+    gradient_op->get_dx(temp,true);
+    diss += 2.0*pow(0.5*(*temp),2);
+    // 6th term: 2e_23^2 = 2*(1/2*(v_z + w_y))^2
+    // v_z
+    find_expansion(grid_type, expan, "v");
+    gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+    gradient_op->get_dz(temp,false);
+    // w_y
+    find_expansion(grid_type, expan, "w");
+    gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+    gradient_op->get_dy(temp,true);
+    diss += 2.0*pow(0.5*(*temp),2);
+    // multiply by 2*mu
+    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) {

src/Science.hpp

@@ -34,6 +34,10 @@ 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);
+// 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
 void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, int Nx, int Ny, int Nz,

src/cases/gravity_current/gravity_current.cpp

@@ -60,7 +60,8 @@ 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?
+bool compute_dissipation;       // Compute dissipation?
+bool compute_BPE;               // Compute background potential energy?
 int iter = 0;                   // Iteration counter
 
 // Maximum squared buoyancy frequency
@@ -188,7 +189,7 @@ class userControl : public BaseCase {
             iter++;
             // Set-up
             if ( iter == 1 ) {
-                if (compute_enstrophy) {
+                if (compute_enstrophy or compute_dissipation) {
                     temp1 = alloc_array(Nx,Ny,Nz);
                 }
                 // Determine last plot if restarting from the dump case
@@ -203,6 +204,7 @@ class userControl : public BaseCase {
             }
 
             /* Calculate and write out useful information */
+
             // Energy (PE assumes density is density anomaly)
             double ke_x, ke_y, ke_z;
             if ( Nx > 1 ) {
@@ -223,6 +225,15 @@ class userControl : public BaseCase {
             if (compute_BPE) {
                 compute_Background_PE(BPE_tot, *tracers[RHO], Nx, Ny, Nz, Lx, Ly, Lz, g, rho_0, iter);
             }
+            // viscous dissipation
+            double diss_tot = 0;
+            double max_diss;
+            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)));
+            }
             // Vorticity / Enstrophy
             double max_vort_x, enst_x_tot;
             double max_vort_y, enst_y_tot;
@@ -277,6 +288,10 @@ class userControl : public BaseCase {
                 if (compute_BPE) {
                     add_diagnostic("BPE_tot", BPE_tot,  header, line);
                 }
+                if (compute_dissipation) {
+                    add_diagnostic("Max_diss", max_diss, header, line);
+                    add_diagnostic("Diss_tot", diss_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);
@@ -427,6 +442,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_dissipation",&compute_dissipation,true,"Calculate dissipation?");
     add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
 
     option_category("Filter options");