Merge pull request #30 from AlexJamesWright/gpu

Gpu updates

Project/CPU/Makefile

@@ -26,7 +26,7 @@ RTFIND_INC_DIR = ./CminpackLibrary/Include
 RTFIND_SRC_DIR = ./CminpackLibrary/Src
 
 # C++ compiler flags
-CXXFLAGS = -Wall -std=c++11 -g -DUSE_MPI=$(USE_MPI) -O3 $(OMP_FLAGS) -Wno-unknown-pragmas
+CXXFLAGS = -Wall -std=c++11 -g -DUSE_MPI=$(USE_MPI) -O0 $(OMP_FLAGS) -Wno-unknown-pragmas
 
 # Sources
 SRCS = main.cc 								\

Project/CPU/Src/interactivePlot.py

@@ -179,19 +179,18 @@ def gatherData(self, states):
                     self.cleanAuxLabels.append(auxLabels[i][:-1])
                 self.cleanAuxLabels.append(auxLabels[-1])
 
-            with suppress(FileNotFoundError):
-                # Grab domain data
-                self.x = np.zeros(c['nx'])
-                self.y = np.zeros(c['ny'])
-                self.z = np.zeros(c['nz'])
-                coords = [self.x, self.y, self.z]
-                print("Fetching domain coordinates...")
-                with open(self.DatDir + 'Domain/domain' + self.appendix +'.dat', 'r') as f:
-                    for coord, (i, line) in zip(coords, enumerate(f)):
-                        temp = line.split()
-                        print(len(temp))
-                        for k, val in enumerate(temp):
-                            coord[k] = float(val)
+#            with suppress(FileNotFoundError):
+#                # Grab domain data
+#                self.x = np.zeros(c['nx'])
+#                self.y = np.zeros(c['ny'])
+#                self.z = np.zeros(c['nz'])
+#                coords = [self.x, self.y, self.z]
+#                print("Fetching domain coordinates...")
+#                with open(self.DatDir + 'Domain/domain' + self.appendix +'.dat', 'r') as f:
+#                    for coord, (i, line) in zip(coords, enumerate(f)):
+#                        temp = line.split()
+#                        for k, val in enumerate(temp):
+#                            coord[k] = float(val)
 
 
 
@@ -259,7 +258,7 @@ def _getYIndexFromLine(self, line, nx, ny):
                 The line number the file pointer is pointing to. We want to know which
                 primitive variable this line's data corresponds to.
             nx: int
-                The total number (incl ghost cells) of domain cells in the x-direction.
+                The total number (incl ghost cells)n of domain cells in the x-direction.
             ny: int
                 The total number (incl ghost cells) of domain cells in the y-direction.
 
@@ -329,10 +328,11 @@ def plotHeatMaps(self, data='prims', color=None, axis=2):
 
             if color==None:
                 color = cm.afmhot
-            surf = ax.imshow(plotVars.T, cmap=color, interpolation='bicubic', aspect='auto')
+            ext = [self.c['xmin'], self.c['xmax'], self.c['ymin'], self.c['ymax']]
+            surf = ax.imshow(plotVars.T, cmap=color, interpolation='bicubic', aspect='auto', origin='lower', extent=ext)
             ax.set_title(r'Time Evolution for {}: $t = {}$'.format(dataLabels[i], c['t']))
-            ax.set_xlim([0, self.c['nx']])
-            ax.set_ylim([0, self.c['ny']])
+            ax.set_xlim([self.c['xmin'], self.c['xmax']])
+            ax.set_ylim([self.c['ymin'], self.c['ymax']])
             ax.set_xlabel(axisLabel1)
             ax.set_ylabel(axisLabel2)
             fig.colorbar(surf, shrink=0.5, aspect=5)
@@ -746,12 +746,4 @@ def plotAdvectionAgainstInitial(self):
 
     Plot = InteractivePlot()
 
-#    Plot.plotSlice()
-#    Plot.plotSingleFluidCurrentSheetAgainstExact()
-#    Plot.plotAdvectionAgainstInitial()
-#    Plot.plotHeatMaps()
-
-    plt.figure()
-    plt.imshow(np.log(Plot.prims[4, :, :, 0].T), extent=[0, 8, 0, 4], origin='lower')
-    plt.show()
-
+#    Plot.plotHeatMaps()

Project/CPU/Src/main.cc

@@ -1,11 +1,12 @@
-// Serial main
-#include "parallelBoundaryConds.h"
+// CPU main
+#include "boundaryConds.h"
 #include "fluxVectorSplitting.h"
-#include "parallelSaveData.h"
+#include "serialSaveData.h"
 #include "simulation.h"
 #include "initFunc.h"
 #include "simData.h"
-#include "RKPlus.h"
+#include "SSP2.h"
+#include "RK2.h"
 #include "Euler.h"
 #include "weno.h"
 
@@ -15,55 +16,50 @@
 using namespace std;
 
 int main(int argc, char *argv[]) {
-
-
+  const double MU(1000);
   // Set up domain
-  int Ng(5);
-  int nx(800);
-  int ny(400);
+  int Ng(4);
+  int nx(64);
+  int ny(16);
   int nz(0);
-  double xmin(0.0);
-  double xmax(8.0);
-  double ymin(0.0);
-  double ymax(4.0);
-  double zmin(0.0);
-  double zmax(1.0);
-  double endTime(30.0);
-  double gamma(2.0);
-  double cfl(0.5);
-  double cp(1);
-  double mu1(-1);
-  double mu2(1);
+  double xmin(-0.5);
+  double xmax(0.5);
+  double ymin(-1.0);
+  double ymax(1.0);
+  double zmin(-1.5);
+  double zmax(1.5);
+  double endTime(0.0005);
+  double cfl(0.1);
+  double gamma(4.0/3.0);
+  double sigma(0);
   bool output(true);
-  int frameSkip(50);
-  int safety(frameSkip);
-  int reportItersPeriod(1);
-  double sigma(50);
-  double nxRanks(4);
+  int safety(180);
+
+  double nxRanks(1);
   double nyRanks(1);
   double nzRanks(1);
 
-  ParallelEnv env(&argc, &argv, nxRanks, nyRanks, nzRanks);
+  SerialEnv env(&argc, &argv, nxRanks, nyRanks, nzRanks);
 
   Data data(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax, endTime, &env,
-            cfl, Ng, gamma, sigma, cp, mu1, mu2, frameSkip, reportItersPeriod);
+            cfl, Ng, gamma, sigma);
 
   // Choose particulars of simulation
-  Euler model(&data);
+  SRMHD model(&data);
 
-  Weno7 weno(&data);
+  Weno3 weno(&data);
 
   FVS fluxMethod(&data, &weno, &model);
 
-  ParallelOutflow bcs(&data, &env);
+  Flow bcs(&data);
 
   Simulation sim(&data, &env);
 
-  FancyMETHODData init(&data);
+  KHInstabilitySingleFluid init(&data, 1);
 
-  RK4 timeInt(&data, &model, &bcs, &fluxMethod);
+  RK2 timeInt(&data, &model, &bcs, &fluxMethod);
 
-  ParallelSaveData save(&data, &env, 0);
+  SerialSaveData save(&data, &env, 0);
 
   // Now objects have been created, set up the simulation
   sim.set(&init, &model, &timeInt, &bcs, &fluxMethod, &save);

Project/CPU/Src/srmhd.cc

@@ -441,6 +441,7 @@ void SRMHD::getPrimitiveVars(double *cons, double *prims, double *aux)
         sol[1] = prims[ID(0, i, j, k)] * aux[ID(0, i, j, k)] /
                  (1 - sol[0]);
 
+
         // Solve residual = 0
         info = __cminpack_func__(hybrd1) (&SRMHDresidual, &args, n, sol, res,
                                           tol, wa, lwa);
@@ -461,8 +462,6 @@ void SRMHD::getPrimitiveVars(double *cons, double *prims, double *aux)
 
 
 
-
-
   // ################################## Smart guessing ########################### //
   // Are there any failures?
   if (fails.size() > 0) {
@@ -507,7 +506,6 @@ void SRMHD::getPrimitiveVars(double *cons, double *prims, double *aux)
     }
   }
 
-
   for (int i(d->is); i < d->ie; i++) {
     for (int j(d->js); j < d->je; j++) {
       for (int k(d->ks); k < d->ke; k++) {

Project/GPU/Include/C2PArgs.h

@@ -21,7 +21,8 @@ class C2PArgs
     bpg,                //!< Blocks per grid
     cellMem,            //!< Memory required for one cell
     Nstreams,           //!< Number of CUDA streams
-    streamWidth;        //!< Number of cells in each stream
+    streamWidth,        //!< Number of cells in each stream
+    nGuessSRMHD;        //!< Number of elements required for the initial guess per cell for the SRMHD model
     double
     //@{
     ** cons_d,

Project/GPU/Include/boundaryConds.h

@@ -21,8 +21,27 @@ class Bcs
         Constructor simply stores the pointer to the Data class.
 
       @param[in] *data pointer to the Data class
+      @param[in] *env pointer to the PlatformEnv class
     */
-    Bcs(Data * data) : data(data) { }
+    Bcs(Data * data, PlatformEnv * env) : data(data)
+    {
+        data->bcsSet = 1;
+    }
+
+    //TODO -- We may not want to allow creation of Bcs object without env in future
+    //! Constructor store data about simulation (needed for domain)
+    /*!
+        Constructor simply stores the pointer to the Data class.
+
+      @param[in] *data pointer to the Data class
+    */
+
+    Bcs(Data * data) : data(data)
+    {
+        data->bcsSet = 1;
+    }
+
+    virtual ~Bcs() { }     //!< Destructor
 
   public:
 
@@ -75,6 +94,8 @@ class Outflow : public Bcs
     */
     Outflow(Data * data) : Bcs(data) { }
 
+    virtual ~Outflow() { }     //!< Destructor
+
     //! Application function
     /*!
         Applies the Outflow boundary conditions to the ghost cells.
@@ -88,6 +109,40 @@ class Outflow : public Bcs
 };
 
 
+//! <b> Out flow boundary conditions for the rotated 2D Brio-Wu </b>
+/*!
+    Using the conventional outflow BCs for the diagonal BW problem results in
+  shocks entering from along the main diagonal. This class deals with these
+  shocks.
+    Using this.apply behaves as if the BW problem has been rotated, as required.
+*/
+class OutflowRotatedBW : public Bcs
+{
+public:
+  //! Constructor
+  /*!
+  Calls constructor of base class to store the pointer to the Data class.
+
+  @param[in] *data pointer to Data class
+  @sa Bcs::Bcs
+  */
+  OutflowRotatedBW(Data * data) : Bcs(data) { }
+
+  virtual ~OutflowRotatedBW() { }     //!< Destructor
+
+  //! Application function
+  /*!
+  Applies the Outflow boundary conditions to the ghost cells.
+
+  @param[in, out] *cons pointer to the conservative (sized) vector
+  @param[in, out] *prims optional pointer to the primitive vector
+  @param[in, out] *aux optional pointer to the primitive vector
+  @sa Bcs::apply
+  */
+  void apply(double * cons, double * prims = NULL, double * aux = NULL);
+};
+
+
 //! <b> Periodic boundary conditions </b>
 /*!
     Flows that exit across one domain boundary re-enter at the opposing
@@ -122,6 +177,8 @@ class Periodic : public Bcs
     */
     Periodic(Data * data) : Bcs(data) { }
 
+    virtual ~Periodic() { }     //!< Destructor
+
     //! Application function
     /*!
         Applies the Periodic boundary conditions to the ghost cells.
@@ -135,16 +192,35 @@ class Periodic : public Bcs
 
 };
 
+//! <b> Flow boundary conditions </b>
 /*!
-  Boundary conditions for the Kelvin Helmholtz instability
-  x-direction is periodic and others are outflow
+    Boundary conditions used for the Kelvin Helmholtz instability. The
+  x-direction is periodic and y- and z-directions are outflow.
 */
 class Flow : public Bcs
 {
 
   public:
+    //! Constructor
+    /*!
+        Calls constructor of base class to store the pointer to the Data class.
+
+      @param[in] *data pointer to Data class
+      @sa Bcs::Bcs
+    */
     Flow(Data * data) : Bcs(data) { }
 
+    virtual ~Flow() { }     //!< Destructor
+
+    //! Application function
+    /*!
+        Applies the Flow boundary conditions to the ghost cells.
+
+      @param[in, out] *cons pointer to the conservative (sized) vector
+      @param[in, out] *prims optional pointer to the primitive vector
+      @param[in, out] *aux optional pointer to the primitive vector
+      @sa Bcs::apply
+    */
     void apply(double * cons, double * prims = NULL, double * aux = NULL);
 
 };