Update to version 0.7

README

@@ -1,4 +1,15 @@
-PENNANT Description:
+The PENNANT Mini-App
+
+Charles R. Ferenbaugh
+Los Alamos National Laboratory
+[email protected]
+
+Version 0.7 -- February 2015
+LA-CC-12-021
+https://github.com/losalamos/PENNANT
+
+
+Description:
 
 PENNANT is an unstructured mesh physics mini-app designed for advanced
 architecture research.  It contains mesh data structures and a few
@@ -11,6 +22,9 @@ PENNANT distribution.
 
 Version Log:
 
+0.7, February 2015:
+     Further optimizations for MPI+OpenMP.
+
 0.6, February 2014:
      First MPI version.  MPI capability is working and mostly
      optimized; MPI+OpenMP is working but needs optimization.

doc/pennantdoc.tex

@@ -19,7 +19,7 @@
 \author{Charles R. Ferenbaugh \\
         Los Alamos National Laboratory \\
         {\tt [email protected]}}
-\date{Version 0.6 -- February 2014 \\
+\date{Version 0.7 -- February 2015 \\
       LA-CC-12-021 \\
       \url{https://github.com/losalamos/PENNANT}}
 \maketitle
@@ -42,7 +42,7 @@ \subsection{Building}
 a simple ``{\tt make}'' command will create a {\tt build} subdirectory and
 build the {\tt pennant} binary in that directory.
 
-PENNANT has been tested under GCC 4.7.2, PGI 13.10, and Intel 13.1.3.
+PENNANT has been tested under GCC 4.7.2, PGI 13.10, and Intel 14.0.4.
 Building under other compilers should require only minor changes.
 
 \subsection{Running tests}
@@ -133,7 +133,7 @@ \subsection{Test inputs and outputs}
 \subsection{Input file parameters}
 
 In most cases, there is no need for users to modify input files.  However,
-here a few parameters that ambitious users might want to know about:
+here are a few parameters that ambitious users might want to know about:
 \begin{description}
     \item[{\tt cstop}]  Stop run when problem reaches given cycle number.
     \item[{\tt tstop}]  Stop run when problem reaches given simulation time.
@@ -209,7 +209,7 @@ \subsection{Mesh data structures}
 PENNANT also stores various mapping arrays from sides to other entity
 types.  These are shown in Figure~\ref{fig:side}.
 Given a side $s$, the following mapping arrays are available:
-\begin{itemize}
+\begin{itemize} \itemsep1pt \parskip0pt
 \item {\tt mapsz} gives the zone $z$ of which $s$ is a subregion.
 \item {\tt mapse} gives the edge $e$ on the boundary of $s$.
 \item {\tt mapsp1} and {\tt mapsp2} give the two mesh points $p_1$ and
@@ -330,7 +330,7 @@ \subsection{Domain decomposition}
     \label{fig:decomp}
 \end{figure}
 
-\begin{enumerate}
+\begin{enumerate} \itemsep1pt \parskip0pt
     \item In {\tt Mesh::parallelGather},
         slave point values are assembled into messages, and sent
         to corresponding {\em proxy} points on the same rank as
@@ -412,7 +412,7 @@ \subsection{Basic hydro algorithms}
 To implement the predictor-corrector scheme, it is necessary to store
 multiple values of some of the problem variables.  This is done
 using the following notation convention:
-\begin{itemize}
+\begin{itemize} \itemsep1pt \parskip0pt
 \item suffix $0$ = the beginning of the timestep (``cycle n'')
 \item suffix $p$ = half-way through the timestep (``cycle n + 1/2'')
 \item no suffix = completion of the timestep    (``cycle n + 1'')
@@ -549,6 +549,9 @@ \section*{Acknowledgements}
 \section{Version Log}
 
 \begin{description}
+\item[0.7] February 2015 \\
+     Further optimizations for MPI+OpenMP.
+
 \item[0.6] February 2014 \\
      First MPI version.  MPI capability is working and mostly
      optimized; MPI+OpenMP is working but needs optimization.

src/Driver.cc

@@ -38,7 +38,7 @@ Driver::Driver(const InputFile* inp, const string& pname)
 
     if (mype == 0) {
         cout << "********************" << endl;
-        cout << "Running PENNANT v0.6" << endl;
+        cout << "Running PENNANT v0.7" << endl;
         cout << "********************" << endl;
         cout << endl;
 
@@ -82,12 +82,13 @@ void Driver::run() {
     time = 0.0;
     cycle = 0;
 
-    double tbegin;
+    double tbegin, tlast;
     if (mype == 0) {
         // get starting timestamp
         struct timeval sbegin;
         gettimeofday(&sbegin, NULL);
         tbegin = sbegin.tv_sec + sbegin.tv_usec * 1.e-6;
+        tlast = tbegin;
     }
 
     // main event loop
@@ -105,12 +106,20 @@ void Driver::run() {
 
         if (mype == 0 &&
                 (cycle == 1 || cycle % dtreport == 0)) {
+            struct timeval scurr;
+            gettimeofday(&scurr, NULL);
+            double tcurr = scurr.tv_sec + scurr.tv_usec * 1.e-6;
+            double tdiff = tcurr - tlast;
+
             cout << scientific << setprecision(5);
             cout << "End cycle " << setw(6) << cycle
                  << ", time = " << setw(11) << time
-                 << ", dt = " << setw(11) << dt << endl;
+                 << ", dt = " << setw(11) << dt
+                 << ", wall = " << setw(11) << tdiff << endl;
             cout << "dt limiter: " << msgdt << endl;
-        }
+
+            tlast = tcurr;
+        } // if mype...
 
     } // while cycle...
 

src/Hydro.cc

@@ -15,6 +15,8 @@
 #include <string>
 #include <vector>
 #include <cmath>
+#include <cstdio>
+#include <cstring>
 #include <algorithm>
 #include <sstream>
 #include <iomanip>
@@ -486,7 +488,7 @@ void Hydro::calcEnergy(
 void Hydro::calcDtCourant(
         const double* zdl,
         double& dtrec,
-        string& msgdtrec,
+        char* msgdtrec,
         const int zfirst,
         const int zlast) {
 
@@ -502,9 +504,7 @@ void Hydro::calcDtCourant(
 
     if (dtnew < dtrec) {
         dtrec = dtnew;
-        ostringstream oss;
-        oss << "Hydro Courant limit for z = " << setw(6) << zmin;
-        msgdtrec = oss.str();
+        snprintf(msgdtrec, 80, "Hydro Courant limit for z = %d", zmin);
     }
 
 }
@@ -515,7 +515,7 @@ void Hydro::calcDtVolume(
         const double* zvol0,
         const double dtlast,
         double& dtrec,
-        string& msgdtrec,
+        char* msgdtrec,
         const int zfirst,
         const int zlast) {
 
@@ -529,9 +529,7 @@ void Hydro::calcDtVolume(
     double dtnew = dtlast * cflv / dvovmax;
     if (dtnew < dtrec) {
         dtrec = dtnew;
-        ostringstream oss;
-        oss << "Hydro dV/V limit for z = " << setw(6) << zmax;
-        msgdtrec = oss.str();
+        snprintf(msgdtrec, 80, "Hydro dV/V limit for z = %d", zmax);
     }
 
 }
@@ -546,16 +544,19 @@ void Hydro::calcDtHydro(
         const int zlast) {
 
     double dtchunk = 1.e99;
-    string msgdtchunk;
+    char msgdtchunk[80];
 
     calcDtCourant(zdl, dtchunk, msgdtchunk, zfirst, zlast);
     calcDtVolume(zvol, zvol0, dtlast, dtchunk, msgdtchunk,
             zfirst, zlast);
-    #pragma omp critical
-    {
-        if (dtchunk < dtrec) {
-            dtrec = dtchunk;
-            msgdtrec = msgdtchunk;
+    if (dtchunk < dtrec) {
+        #pragma omp critical
+        {
+            // redundant test needed to avoid race condition
+            if (dtchunk < dtrec) {
+                dtrec = dtchunk;
+                strncpy(msgdtrec, msgdtchunk, 80);
+            }
         }
     }
 
@@ -568,7 +569,7 @@ void Hydro::getDtHydro(
 
     if (dtrec < dtnew) {
         dtnew = dtrec;
-        msgdtnew = msgdtrec;
+        msgdtnew = string(msgdtrec);
     }
 
 }
@@ -577,6 +578,6 @@ void Hydro::getDtHydro(
 void Hydro::resetDtHydro() {
 
     dtrec = 1.e99;
-    msgdtrec = "Hydro default";
+    strcpy(msgdtrec, "Hydro default");
 
 }

src/Hydro.hh

@@ -50,7 +50,7 @@ public:
     std::vector<double> bcy;    // y values of y-plane fixed boundaries
 
     double dtrec;               // maximum timestep for hydro
-    std::string msgdtrec;       // message:  reason for dtrec
+    char msgdtrec[80];          // message:  reason for dtrec
 
     double2* pu;       // point velocity
     double2* pu0;      // point velocity, start of cycle
@@ -168,7 +168,7 @@ public:
     void calcDtCourant(
             const double* zdl,
             double& dtrec,
-            std::string& msgdtrec,
+            char* msgdtrec,
             const int zfirst,
             const int zlast);
 
@@ -177,7 +177,7 @@ public:
             const double* zvol0,
             const double dtlast,
             double& dtrec,
-            std::string& msgdtrec,
+            char* msgdtrec,
             const int zfirst,
             const int zlast);
 

src/Memory.hh

@@ -14,6 +14,9 @@
 #define MEMORY_HH_
 
 #include <cstdlib>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 
 // Namespace Memory provides functions to allocate and free memory.
@@ -25,12 +28,20 @@ namespace Memory {
 
 template<typename T>
 inline T* alloc(const int count) {
+#if defined(_OPENMP) && defined(__INTEL_COMPILER)
+    return (T*) kmp_malloc(count * sizeof(T));
+#else
     return (T*) std::malloc(count * sizeof(T));
+#endif
 }
 
 template<typename T>
 inline void free(T* ptr) {
+#if defined(_OPENMP) && defined(__INTEL_COMPILER)
+    kmp_free(ptr);
+#else
     std::free(ptr);
+#endif
 }
 
 };  // namespace Memory

src/Mesh.cc

@@ -159,9 +159,9 @@ void Mesh::init() {
 
 
 void Mesh::initSides(
-        std::vector<int>& cellstart,
-        std::vector<int>& cellsize,
-        std::vector<int>& cellnodes) {
+        const vector<int>& cellstart,
+        const vector<int>& cellsize,
+        const vector<int>& cellnodes) {
 
     mapsp1 = Memory::alloc<int>(nums);
     mapsp2 = Memory::alloc<int>(nums);
@@ -285,12 +285,12 @@ void Mesh::initInvMap() {
 
 
 void Mesh::initParallel(
-        std::vector<int>& slavemstrpes,
-        std::vector<int>& slavemstrcounts,
-        std::vector<int>& slavepoints,
-        std::vector<int>& masterslvpes,
-        std::vector<int>& masterslvcounts,
-        std::vector<int>& masterpoints) {
+        const vector<int>& slavemstrpes,
+        const vector<int>& slavemstrcounts,
+        const vector<int>& slavepoints,
+        const vector<int>& masterslvpes,
+        const vector<int>& masterslvcounts,
+        const vector<int>& masterpoints) {
     if (Parallel::numpe == 1) return;
 
     nummstrpe = slavemstrpes.size();
@@ -492,8 +492,10 @@ void Mesh::calcVols(
 
     } // for s
 
-    #pragma omp atomic
-    numsbad += count;
+    if (count > 0) {
+        #pragma omp atomic
+        numsbad += count;
+    }
 
 }
 

src/Mesh.hh

@@ -115,9 +115,9 @@ public:
 
     // populate mapping arrays
     void initSides(
-            std::vector<int>& cellstart,
-            std::vector<int>& cellsize,
-            std::vector<int>& cellnodes);
+            const std::vector<int>& cellstart,
+            const std::vector<int>& cellsize,
+            const std::vector<int>& cellnodes);
     void initEdges();
 
     // populate chunk information
@@ -127,12 +127,12 @@ public:
     void initInvMap();
 
     void initParallel(
-            std::vector<int>& slavemstrpes,
-            std::vector<int>& slavemstrcounts,
-            std::vector<int>& slavepoints,
-            std::vector<int>& masterslvpes,
-            std::vector<int>& masterslvcounts,
-            std::vector<int>& masterpoints);
+            const std::vector<int>& slavemstrpes,
+            const std::vector<int>& slavemstrcounts,
+            const std::vector<int>& slavepoints,
+            const std::vector<int>& masterslvpes,
+            const std::vector<int>& masterslvcounts,
+            const std::vector<int>& masterpoints);
 
     // write mesh statistics
     void writeStats();