Initial commits for instrumentation/scaling benchmarking

src/cases/bench_dipole.cpp

@@ -12,6 +12,9 @@
 #include <random/normal.h>
 #include <vector>
 #include "../Science.hpp"
+#include "../multigrid.hpp"
+#include "../grad.hpp"
+#include "../gmres.hpp"
 
 using namespace std;
 using namespace TArrayn;
@@ -310,6 +313,14 @@ int main(int argc, char ** argv) {
    ppois.do_run(0.001);
    now = MPI_Wtime();
    if (master()) fprintf(stderr,"Total runtime (%d iterations) complete in %gs (%gs per)\n",mycode.itercount,now-start_time,(now-start_time)/mycode.itercount);
+   if (master()) fprintf(stderr,"  %d invocations of coarse solve, for %gs (%gs per)\n",
+                  coarse_solve_count,coarse_solve_time,coarse_solve_time/coarse_solve_count);
+   if (master()) fprintf(stderr,"  %d red black smoothings, for %gs (%gs per)\n",
+                  redblack_count, redblack_time, redblack_time/redblack_count);
+   if (master()) fprintf(stderr,"  %d FD operator applications, for %gs (%gs per)\n",
+                  apply_op_count, apply_op_time, apply_op_time/apply_op_count);
+   if (master()) fprintf(stderr,"  %gs spent in spectral derivatives\n",deriv_time);
+   if (master()) fprintf(stderr,"  %gs spent in gmres-lapack internals\n",gmres_lapack_time);
    MPI_Finalize();
    return 0;
 }

src/cases/tankrho/spins.conf

+restart = false
+Lx = 1.0
+Ly = 1.0
+Lz = 1.0
+Nx = 1024 
+Ny = 1
+Nz = 1024
+
+rho_0 = 1.0
+delta_rho = 0.01
+pyc_asym = 0.0
+pyc_sep_perc = 0.0
+h_perc =0.1
+h_mix_perc = 0.1
+delta_x = 0.01
+Lmix = 0.1
+Hmix = 0.1
+
+hill_height = 0.1
+hill_centre = 0.5
+hill_width  = 0.1
+
+visco = 0.0001
+diffu_rho = 0.00001
+
+plot_interval = 1.0
+final_time = 0.5
+
+perturb  = 0.0001
+

src/gmres.cpp

 #include "gmres.hpp"
 
+double gmres_lapack_time = 0;
 // Defines whether GMRES should be chatty about output, for debugging purposes

src/gmres.hpp

@@ -4,6 +4,7 @@
 #include <blitz/array.h>
 #include <stdio.h>
 #include <iostream>
+#include <mpi.h>
 //#include <mkl_lapack.h>
 /* LAPACK function prototype */
 extern "C" {
@@ -13,6 +14,8 @@ extern "C" {
                int * LDB, double * S, double * RCOND, int * RANK, double * WORK,
                int * LWORK, int * IWORK, int * INFO);
 }
+
+extern double gmres_lapack_time;
 /* gmres.hpp -- header class for gmres template class, which abstracts out the matrix-vector
    multiply, preconditioning, and dot products of GMRES to a user-specified type conforming
    to a provided interface.
@@ -313,9 +316,11 @@ template <class Controller> class GMRES_Solver {
 //               dgels_("N",&M,&N,&NRHS,hess_copy.data(),&LDA,
 //                     rhs_vec.data(), &LDB, lapack_workspace, &lwork_size,
 //                     &INFO);
+               double now = MPI_Wtime();
                dgelsd_(&M,&N,&NRHS,hess_copy.data(),&LDA,rhs_vec.data(),
                      &LDB, svd_vec.data(), &RCOND, &RANK, lapack_workspace,
                     &lwork_size,lapack_iworkspace,&INFO); 
+               gmres_lapack_time += (MPI_Wtime() - now);
 //               fprintf(stderr,"%s:%d RANK %d\n",__FILE__,__LINE__,RANK); fflush(stderr);
                /* rhs_vec.data() contains the answer and remainder */
 //               std::cerr << svd_vec;

src/grad.cpp

@@ -5,6 +5,8 @@
 #include "T_util.hpp"
 #include <stdio.h>
 
+double deriv_time = 0;
+
 using namespace Transformer;
 namespace TArrayn {
 
@@ -135,6 +137,7 @@ namespace TArrayn {
    /* Private function to save boilerplate of calculating d/d[abc] with the
       right kind of transform */
    void Grad::calculate_d(S_EXP type, Trans1D * tform, blitz::Array<double,3> * restrict dest) {
+      double now = MPI_Wtime();
       switch(type) {
          case COSINE:
             deriv_dct(*my_array,*tform,*dest);
@@ -152,6 +155,7 @@ namespace TArrayn {
             fprintf(stderr,"Invalid expansion type (%d) in grad!\n",type);
             abort();
       }
+      deriv_time += (MPI_Wtime() - now);
    }
 
    void Grad::do_diff(DTArray * restrict dx_out, Dimension inDim) {

src/grad.hpp

@@ -14,6 +14,8 @@
 #include "Parformer.hpp"
 #include <blitz/array.h>
 
+extern double deriv_time;
+
 namespace TArrayn {
    /* Include in the TArrayn namespace, since this is used much like
       deriv_cheb / deriv_fft / etc */

src/multigrid.cpp

@@ -10,6 +10,16 @@
 using namespace blitz;
 using namespace std;
 
+double coarse_solve_time = 0;
+int coarse_solve_count = 0;
+
+double redblack_time = 0;
+int redblack_count = 0;
+
+double apply_op_time = 0;
+int apply_op_count = 0;
+
+
 //#define COARSE_GRID_SIZE 512
 #define COARSE_GRID_SIZE 8
 #define SYNC(__x__) { for (int _i = 0; _i < nproc; _i++) { \
@@ -1830,10 +1840,13 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
             u_ptr = u.data();
             f_ptr = f.data();
          }
+         double now = MPI_Wtime();
          int retval = umfpack_di_solve(sys,A_cols,A_rows,A_double,
                               u_ptr, f_ptr,
                               numeric_factor,
                               Control,Info);
+         coarse_solve_time += (MPI_Wtime() - now);
+         coarse_solve_count += 1;
          assert(!retval);
          if (indefinite_problem) {
 //            fprintf(stderr,"CG: True u\n");
@@ -1868,10 +1881,16 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
    /* Perform one red-black relaxtion pass */
 //   if (master()) fprintf(stderr,"Pre-cyle one\n");
    if (pre >= 1) {
+      double now = MPI_Wtime();
       do_redblack(f,correction);
+      double later = MPI_Wtime();
+      redblack_time += (later - now);
+      redblack_count++;
    //cout << "New correction: " << correction;
    /* And apply the operator to find the defect */
       apply_operator(correction,defect);
+      apply_op_count++;
+      apply_op_time += (MPI_Wtime() - later);
       defect = -(defect - f);
       u = correction;
    /* Perform a second relaxation pass */
@@ -1887,11 +1906,17 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
    }
    for (int i = 1; i < pre; i++) {
 //      if (master()) fprintf(stderr,"Doing RB precycle %d\n",i);
+      double now = MPI_Wtime();
       do_redblack(defect,correction);
+      double later = MPI_Wtime();
+      redblack_time += (later - now);
+      redblack_count++;
 //      if (master()) fprintf(stderr,"RB precycle %d complete\n",i);
 //      if (master()) fprintf(stderr,"correction:\n");
 //      SYNC(cerr << correction; cerr.flush());
       apply_operator(correction,f);
+      apply_op_count++;
+      apply_op_time += (MPI_Wtime() - later);
 //      if (master()) fprintf(stderr,"f:\n");
 //      SYNC(cerr << f; cerr.flush());
       u = u+correction;
@@ -1936,7 +1961,12 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
          rebalance_array(coarse_u,local_coarse_2d,my_comm,UCoarse);
          interpolate_grid(correction);
 //         cout << "CG: New correction: " << correction;
-         apply_operator(correction,f);
+         {
+            double now = MPI_Wtime();
+            apply_operator(correction,f);
+            apply_op_count++;
+            apply_op_time += (MPI_Wtime() - now);
+         }
          defect = defect-f-coarse_extra_out;
          u = u+correction;
          extra_out = coarse_extra_out;
@@ -1965,7 +1995,12 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
          rebalance_array(coarse_u,local_coarse_2d,my_comm,UCoarse);
          interpolate_grid(correction);
 //         cout << "CG: New correction: " << correction;
-         apply_operator(correction,f);
+         {
+            double now = MPI_Wtime();
+            apply_operator(correction,f);
+            apply_op_count++;
+            apply_op_time += (MPI_Wtime() - now);
+         }
 //         SYNC(cerr << f; cerr.flush()); MPI_Finalize(); exit(1);
          defect = defect-f-coarse_extra_out;
          u = u+correction;
@@ -1974,9 +2009,15 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
          /* Now, two more relaxation steps */
 //         if (master()) fprintf(stderr,"F[%d] mid-cycle corrections\n",size_x);
          for (int i = 0; i < mid; i++) {
+            double now = MPI_Wtime();
             do_redblack(defect,correction);
       //      cout << "New correction: " << correction;
+            double later = MPI_Wtime();
+            redblack_count++;
+            redblack_time += (later - now);
             apply_operator(correction,f);
+            apply_op_count++;
+            apply_op_time += (MPI_Wtime() - later);
             defect = defect-f;
             u = u+correction;
          }
@@ -1998,7 +2039,12 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
          MPI_Bcast(&coarse_extra_out,1,MPI_DOUBLE,0,my_comm);
          rebalance_array(coarse_u,local_coarse_2d,my_comm,UCoarse);
          interpolate_grid(correction);
-         apply_operator(correction,f);
+         {
+            double now = MPI_Wtime();
+            apply_operator(correction,f);
+            apply_op_count++;
+            apply_op_time += (MPI_Wtime() - now);
+         }
          defect = defect-f-coarse_extra_out;
          u = u+correction;
 //         fprintf(stderr,"F(2)[%d] average: %g vs %g\n",size_x,pvsum(u,my_comm)/(size_x*size_z),extra_in);
@@ -2016,9 +2062,15 @@ void MG_Solver::_cycle(CYCLE_TYPE cycle, Array<double,2> & f, Array<double,2> &
 
    /* Now, two more relaxation steps */
    for (int i = 0; i < post; i++) {
+      double now = MPI_Wtime();
       do_redblack(defect,correction);
+      double later = MPI_Wtime();
+      redblack_count++;
+      redblack_time += (later - now);
 //      cout << "New correction: " << correction;
       apply_operator(correction,f);
+      apply_op_count++;
+      apply_op_time += (MPI_Wtime() - later);
       defect = defect-f;
       u = u+correction;
    }

src/multigrid.hpp

@@ -8,6 +8,10 @@
 /* Header file for our required 2D multigrid solver, operating on a (split) grid
    with the numerical Jacobian entries (from the spectral model) providing varying
    coefficients */
+
+extern double coarse_solve_time, redblack_time, apply_op_time;
+extern int coarse_solve_count, redblack_count, apply_op_count;
+
 enum SYM_TYPE {
    SYM_NONE=0, // No symmetry
    SYM_PERIODIC, // Periodic