First part of the flux computation in 3D (not finished!)

src/GksGpu/CellUpdate/CellUpdate.cu

@@ -22,7 +22,7 @@ __host__ __device__ inline void cellUpdateFunction( DataBaseStruct dataBase, Par
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-void CellUpdate::updateCells( SPtr<DataBase> dataBase, Parameters parameters, uint level )
+void CellUpdate::run( SPtr<DataBase> dataBase, Parameters parameters, uint level )
 {
     CudaUtility::CudaGrid grid( dataBase->perLevelCount[ level ].numberOfBulkCells, 32 );
 
@@ -33,7 +33,7 @@ void CellUpdate::updateCells( SPtr<DataBase> dataBase, Parameters parameters, ui
                parameters,
                dataBase->perLevelCount[ level ].startOfCells );
 
-    getLastCudaError("CellUpdate::updateCells( SPtr<DataBase> dataBase, Parameters parameters, uint level )");
+    getLastCudaError("CellUpdate::run( SPtr<DataBase> dataBase, Parameters parameters, uint level )");
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

src/GksGpu/CellUpdate/CellUpdate.h

@@ -13,9 +13,9 @@ class VF_PUBLIC CellUpdate
 {
 public:
 
-    static void updateCells( SPtr<DataBase> dataBase, 
-                             Parameters parameters, 
-                             uint level );
+    static void run( SPtr<DataBase> dataBase, 
+                     Parameters parameters, 
+                     uint level );
 };
 
 #endif

src/GksGpu/FluxComputation/FluxComputation.cu

+#include "FluxComputation.h"
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include "Core/PointerDefinitions.h"
+#include "Core/RealConstants.h"
+
+#include "DataBase/DataBaseStruct.h"
+
+#include "Definitions/PassiveScalar.h"
+
+#include "FlowStateData/FlowStateData.cuh"
+
+#include "CudaUtility/CudaRunKernel.hpp"
+
+__global__                 void fluxKernel  ( DataBaseStruct dataBase, Parameters parameters, char direction, uint startIndex, uint numberOfEntities );
+
+__host__ __device__ inline void fluxFunction( DataBaseStruct dataBase, Parameters parameters, char direction, uint startIndex, uint index );
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void FluxComputation::run( SPtr<DataBase> dataBase, Parameters parameters, uint level )
+{
+    {
+        CudaUtility::CudaGrid grid(dataBase->perLevelCount[level].numberOfFacesX, 32);
+
+        runKernel(fluxKernel,
+            fluxFunction,
+            dataBase->getDeviceType(), grid,
+            dataBase->toStruct(),
+            parameters,
+            'x',
+            dataBase->perLevelCount[level].startOfFacesX);
+
+        getLastCudaError("FluxComputation::run( SPtr<DataBase> dataBase, Parameters parameters, 'x', uint level )");
+    }
+    {
+        CudaUtility::CudaGrid grid(dataBase->perLevelCount[level].numberOfFacesY, 32);
+
+        runKernel(fluxKernel,
+            fluxFunction,
+            dataBase->getDeviceType(), grid,
+            dataBase->toStruct(),
+            parameters,
+            'y',
+            dataBase->perLevelCount[level].startOfFacesY);
+
+        getLastCudaError("FluxComputation::run( SPtr<DataBase> dataBase, Parameters parameters, 'y', uint level )");
+    }
+    {
+        CudaUtility::CudaGrid grid(dataBase->perLevelCount[level].numberOfFacesZ, 32);
+
+        runKernel(fluxKernel,
+            fluxFunction,
+            dataBase->getDeviceType(), grid,
+            dataBase->toStruct(),
+            parameters,
+            'z',
+            dataBase->perLevelCount[level].startOfFacesZ);
+
+        getLastCudaError("FluxComputation::run( SPtr<DataBase> dataBase, Parameters parameters, 'z', uint level )");
+    }
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+__global__ void fluxKernel(DataBaseStruct dataBase, Parameters parameters, char direction, uint startIndex, uint numberOfEntities)
+{
+    uint index = blockIdx.x * blockDim.x + threadIdx.x;
+
+    if( index > numberOfEntities ) return;
+
+    fluxFunction( dataBase, parameters, direction, startIndex, index );
+}
+
+__host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters parameters, char direction, uint startIndex, uint index)
+{
+    uint faceIndex = startIndex + index;
+
+    //////////////////////////////////////////////////////////////////////////
+
+
+
+}

src/GksGpu/FluxComputation/FluxComputation.h

+#ifndef  FluxComputation_H
+#define  FluxComputation_H
+
+#include "VirtualFluidsDefinitions.h"
+
+#include "Core/PointerDefinitions.h"
+#include "Core/DataTypes.h"
+
+#include "DataBase/DataBase.h"
+#include "Parameters/Parameters.h"
+
+class VF_PUBLIC FluxComputation
+{
+public:
+
+    static void run( SPtr<DataBase> dataBase, 
+                     Parameters parameters, 
+                     uint level );
+};
+
+#endif

src/GksGpu/FluxComputation/Reconstruction.cuh

+#ifndef Reconstruction_CUH
+#define Reconstruction_CUH
+
+#include "VirtualFluidsDefinitions.h"
+
+#include "Core/DataTypes.h"
+
+#include "DataBase/DataBase.h"
+#include "Parameters/Parameters.h"
+
+#include "FlowStateData/FlowStateData.cuh"
+#include "FlowStateData/AccessDeviceData.cuh"
+
+__host__ __device__ inline void getCellIndicesN ( const uint faceIndex,
+                                                  const DataBaseStruct& dataBase,
+                                                  uint& posCellIndexN,
+                                                  uint& negCellIndexN )
+{
+    posCellIndexN = dataBase.faceToCell[ POS_CELL( faceIndex, dataBase.numberOfFaces ) ];
+    negCellIndexN = dataBase.faceToCell[ NEG_CELL( faceIndex, dataBase.numberOfFaces ) ];
+}
+
+__host__ __device__ inline void getCellIndicesTX( const uint faceIndex,
+                                                  const DataBaseStruct& dataBase,
+                                                  const uint posCellIndexN,
+                                                  const uint negCellIndexN,
+                                                  uint* posCellIndexTX,
+                                                  uint* negCellIndexTX )
+{
+    posCellIndexTX[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 0, dataBase.numberOfCells ) ];
+    posCellIndexTX[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 0, dataBase.numberOfCells ) ];
+
+    negCellIndexTX[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 1, dataBase.numberOfCells ) ];
+    negCellIndexTX[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 1, dataBase.numberOfCells ) ];
+}
+
+__host__ __device__ inline void getCellIndicesTY( const uint faceIndex,
+                                                  const DataBaseStruct& dataBase,
+                                                  const uint posCellIndexN,
+                                                  const uint negCellIndexN,
+                                                  uint* posCellIndexTY,
+                                                  uint* negCellIndexTY )
+{
+    posCellIndexTY[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 2, dataBase.numberOfCells ) ];
+    posCellIndexTY[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 2, dataBase.numberOfCells ) ];
+
+    negCellIndexTY[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 3, dataBase.numberOfCells ) ];
+    negCellIndexTY[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 3, dataBase.numberOfCells ) ];
+}
+
+__host__ __device__ inline void getCellIndicesTZ( const uint faceIndex,
+                                                  const DataBaseStruct& dataBase,
+                                                  const uint posCellIndexN,
+                                                  const uint negCellIndexN,
+                                                  uint* posCellIndexTZ,
+                                                  uint* negCellIndexTZ )
+{
+    posCellIndexTZ[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 4, dataBase.numberOfCells ) ];
+    posCellIndexTZ[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 4, dataBase.numberOfCells ) ];
+
+    negCellIndexTZ[0] = dataBase.cellToCell[ CELL_TO_CELL( posCellIndexN, 5, dataBase.numberOfCells ) ];
+    negCellIndexTZ[1] = dataBase.cellToCell[ CELL_TO_CELL( negCellIndexN, 5, dataBase.numberOfCells ) ];
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+__host__ __device__ inline void computeFaceCons( const ConservedVariables& posCons,
+                                                 const ConservedVariables& negCons,
+                                                 ConservedVariables& faceCons )
+{
+    faceCons.rho  = c1o2 * ( negCons.rho  + posCons.rho  );
+    faceCons.rhoU = c1o2 * ( negCons.rhoU + posCons.rhoU );
+    faceCons.rhoV = c1o2 * ( negCons.rhoV + posCons.rhoV );
+    faceCons.rhoW = c1o2 * ( negCons.rhoW + posCons.rhoW );
+    faceCons.rhoE = c1o2 * ( negCons.rhoE + posCons.rhoE );
+#ifdef USE_PASSIVE_SCALAR
+	faceCons.rhoS = c1o2 * ( negCons.rhoS + posCons.rhoS );
+#endif // USE_PASSIVE_SCALAR
+}
+
+__host__ __device__ inline void computeGradN( const Parameters& parameters,
+                                              const ConservedVariables& posCons,
+                                              const ConservedVariables& negCons,
+                                              const PrimitiveVariables& facePrim,
+                                              ConservedVariables& gradN )
+{
+    gradN.rho  = ( posCons.rho  - negCons.rho  ) / ( parameters.dx * facePrim.rho );
+    gradN.rhoU = ( posCons.rhoU - negCons.rhoU ) / ( parameters.dx * facePrim.rho );
+    gradN.rhoV = ( posCons.rhoV - negCons.rhoV ) / ( parameters.dx * facePrim.rho );
+    gradN.rhoW = ( posCons.rhoW - negCons.rhoW ) / ( parameters.dx * facePrim.rho );
+    gradN.rhoE = ( posCons.rhoE - negCons.rhoE ) / ( parameters.dx * facePrim.rho );
+#ifdef USE_PASSIVE_SCALAR
+	gradN.rhoS = ( posCons.rhoS - negCons.rhoS ) / ( parameters.dx * facePrim.rho );
+#endif // USE_PASSIVE_SCALAR
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+__host__ __device__ inline void reconstructFiniteDifferences( const uint faceIndex,
+                                                              const DataBaseStruct& dataBase,
+                                                              const Parameters& parameters,
+                                                              const char direction,
+                                                              ConservedVariables& gradN,
+                                                              ConservedVariables& gradT1,
+                                                              ConservedVariables& gradT2,
+                                                              PrimitiveVariables& facePrim )
+{
+    uint posCellIndexN, negCellIndexN;
+
+    getCellIndicesN( faceIndex, dataBase, posCellIndexN, negCellIndexN );
+    
+    {
+
+        ConservedVariables posCons, negCons, faceCons;
+
+        readCellData(posCellIndexN, dataBase, posCons);
+        readCellData(negCellIndexN, dataBase, negCons);
+        
+        computeFaceCons(posCons, negCons, faceCons);
+
+        facePrim = toPrimitiveVariables( faceCons, parameters.K );
+
+        computeGradN( parameters, posCons, negCons, facePrim, gradN );
+    }
+
+    {
+        uint posCellIndexT1[2];
+        uint negCellIndexT1[2];
+    
+        if( direction == 'x' ) getCellIndicesTY(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT1, negCellIndexT1);
+        if( direction == 'y' ) getCellIndicesTZ(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT1, negCellIndexT1);
+        if( direction == 'z' ) getCellIndicesTX(faceIndex, dataBase, posCellIndexN, negCellIndexN, posCellIndexT1, negCellIndexT1);
+
+        ConservedVariables cons;
+
+        //////////////////////////////////////////////////////////////////////////
+        {
+            readCellData(posCellIndexT1[0], dataBase, cons);
+
+            gradN.rho  += c1o2 * cons.rho;
+            gradN.rhoU += c1o2 * cons.rhoU;
+            gradN.rhoV += c1o2 * cons.rhoV;
+            gradN.rhoW += c1o2 * cons.rhoW;
+            gradN.rhoE += c1o2 * cons.rhoE;
+        #ifdef USE_PASSIVE_SCALAR
+            gradN.rhoS += c1o2 * cons.rhoS;
+        #endif // USE_PASSIVE_SCALAR
+        }
+        {
+            readCellData(posCellIndexT1[1], dataBase, cons);
+
+            gradN.rho  += c1o2 * cons.rho;
+            gradN.rhoU += c1o2 * cons.rhoU;
+            gradN.rhoV += c1o2 * cons.rhoV;
+            gradN.rhoW += c1o2 * cons.rhoW;
+            gradN.rhoE += c1o2 * cons.rhoE;
+        #ifdef USE_PASSIVE_SCALAR
+            gradN.rhoS += c1o2 * cons.rhoS;
+        #endif // USE_PASSIVE_SCALAR
+        }
+        //////////////////////////////////////////////////////////////////////////
+        {
+            readCellData(negCellIndexT1[0], dataBase, cons);
+
+            gradN.rho  -= c1o2 * cons.rho;
+            gradN.rhoU -= c1o2 * cons.rhoU;
+            gradN.rhoV -= c1o2 * cons.rhoV;
+            gradN.rhoW -= c1o2 * cons.rhoW;
+            gradN.rhoE -= c1o2 * cons.rhoE;
+        #ifdef USE_PASSIVE_SCALAR
+            gradN.rhoS -= c1o2 * cons.rhoS;
+        #endif // USE_PASSIVE_SCALAR
+        }
+        {
+            readCellData(negCellIndexT1[1], dataBase, cons);
+
+            gradN.rho  -= c1o2 * cons.rho;
+            gradN.rhoU -= c1o2 * cons.rhoU;
+            gradN.rhoV -= c1o2 * cons.rhoV;
+            gradN.rhoW -= c1o2 * cons.rhoW;
+            gradN.rhoE -= c1o2 * cons.rhoE;
+        #ifdef USE_PASSIVE_SCALAR
+            gradN.rhoS -= c1o2 * cons.rhoS;
+        #endif // USE_PASSIVE_SCALAR
+        }
+        //////////////////////////////////////////////////////////////////////////
+        {
+            gradN.rho  /= two * parameters.dx * facePrim.rho;
+            gradN.rhoU /= two * parameters.dx * facePrim.rho;
+            gradN.rhoV /= two * parameters.dx * facePrim.rho;
+            gradN.rhoW /= two * parameters.dx * facePrim.rho;
+            gradN.rhoE /= two * parameters.dx * facePrim.rho;
+        #ifdef USE_PASSIVE_SCALAR
+            gradN.rhoS /= two * parameters.dx * facePrim.rho;
+        #endif // USE_PASSIVE_SCALAR
+        }
+    }
+}
+
+
+
+
+
+
+#endif
\ No newline at end of file

src/GksGpu/FluxEvaluation/utilityKernels.cu

-#include <cstdio>
-#include <cuda_runtime.h>
-#include <device_launch_parameters.h>
-
-#include "DataBase/DataBaseStruct.h"
-
-#define THREADS_PER_BLOCK 128
-
-__global__ void dummyKernel();
-
-__global__ void debugKernel( const DataBaseStruct dataBase );
-
-//////////////////////////////////////////////////////////////////////////
-
-__global__ void dummyKernel(  )
-{
-    printf("I am thread %d.\n", threadIdx.x);
-}
-
-__global__ void debugKernel( const DataBaseStruct dataBase )
-{   
-    if( threadIdx.x + blockIdx.x == 0 ){
-        printf("numberOfCells     : %d\n", dataBase.numberOfCells     );
-        printf("numberOfFaces     : %d\n", dataBase.numberOfFaces     );
-        printf("\n");
-        printf("\n");
-        printf("faceToCell: %p\n", dataBase.faceToCell);
-        printf("faceCenter: %p\n", dataBase.faceCenter);
-        printf("cellCenter: %p\n", dataBase.cellCenter);
-        printf("\n");
-        printf("data      : %p\n", dataBase.data      );
-        printf("dataNew   : %p\n", dataBase.dataUpdate);
-    }
-}

src/GksGpu/TimeStepping/NestedTimeStep.cpp

@@ -31,6 +31,6 @@ void TimeStepping::nestedTimeStep( SPtr<DataBase> dataBase,
 
     //runFluxKernel( dataBase, parameters, type, level ); getLastCudaError();
 
-    CellUpdate::updateCells( dataBase, parameters, level );
+    CellUpdate::run( dataBase, parameters, level );
 }