diff --git a/src/GksGpu/BoundaryConditions/Open.cu b/src/GksGpu/BoundaryConditions/Open.cu
index cbd803a5f0844dadeb416b4ce8c9ebcc4840417b..e35ff0276d30225de3768e31b7ddc5ecc30bfe08 100644
--- a/src/GksGpu/BoundaryConditions/Open.cu
+++ b/src/GksGpu/BoundaryConditions/Open.cu
@@ -114,11 +114,22 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
real p1 = c1o2 * domainCellPrim.rho / domainCellPrim.lambda;
real p2 = c1o2 * secondCellPrim.rho / secondCellPrim.lambda;
+ real p0 = c1o2 * boundaryCondition.prim.rho / boundaryCondition.prim.lambda;
+
//////////////////////////////////////////////////////////////////////////
- if( sign > zero )
+ //if( sign > zero )
//if( p2 > p1 )
+ if( p1 > p0 )
+ {
ghostCellData = domainCellData;
+ //ghostCellData = two * domainCellData + ( - one ) * secondCellData;
+
+
+ ghostCellData.rhoU = zero;
+ ghostCellData.rhoV = zero;
+ ghostCellData.rhoW = zero;
+ }
else
{
PrimitiveVariables ghostCellPrim = boundaryCondition.prim;
@@ -127,6 +138,14 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
ghostCellPrim.V = domainCellPrim.V;
ghostCellPrim.W = domainCellPrim.W;
+ //ghostCellPrim.U = p0/p1;
+ //ghostCellPrim.V = p0/p1;
+ //ghostCellPrim.W = p0/p1;
+
+ //ghostCellPrim.U = two * domainCellPrim.U - secondCellPrim.U;
+ //ghostCellPrim.V = two * domainCellPrim.V - secondCellPrim.V;
+ //ghostCellPrim.W = two * domainCellPrim.W - secondCellPrim.W;
+
real velocity = sqrt( ghostCellPrim.U * ghostCellPrim.U + ghostCellPrim.V * ghostCellPrim.V + ghostCellPrim.W * ghostCellPrim.W );
if( velocity > ten )
diff --git a/src/GksGpu/BoundaryConditions/Pressure.cu b/src/GksGpu/BoundaryConditions/Pressure.cu
index c6848c85676a16b0b0ad8b170a53754fc5491907..2c028a7d6e13a381a2a88868c6673de41ba2f4b6 100644
--- a/src/GksGpu/BoundaryConditions/Pressure.cu
+++ b/src/GksGpu/BoundaryConditions/Pressure.cu
@@ -103,6 +103,7 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
}
}
+ //ghostCellPrim.rho = two * domainCellPrim.rho - secondCellPrim.rho;
ghostCellPrim.U = two * domainCellPrim.U - secondCellPrim.U;
ghostCellPrim.V = two * domainCellPrim.V - secondCellPrim.V;
ghostCellPrim.W = two * domainCellPrim.W - secondCellPrim.W;
@@ -112,18 +113,12 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
ghostCellPrim.S_2 = two * domainCellPrim.S_2 - secondCellPrim.S_2;
#endif // USE_PASSIVE_SCALAR
- // ghostCellPrim.U = domainCellPrim.U ;
- // ghostCellPrim.V = domainCellPrim.V ;
- // ghostCellPrim.W = domainCellPrim.W ;
- // ghostCellPrim.lambda = domainCellPrim.lambda;
- //#ifdef USE_PASSIVE_SCALAR
- // ghostCellPrim.S = domainCellPrim.S ;
- //#endif // USE_PASSIVE_SCALAR
-
real rho0 = ( two * boundaryCondition.p0 * c1o2 * ( domainCellPrim.lambda + ghostCellPrim.lambda ) );
-
ghostCellPrim.rho = two * rho0 - domainCellPrim.rho;
+
+ //real lambda0 = ( c1o2 * ( domainCellPrim.rho + ghostCellPrim.rho ) * c1o2 / boundaryCondition.p0 );
+ //ghostCellPrim.lambda = two * lambda0 - domainCellPrim.lambda;
}
{
diff --git a/src/GksGpu/BoundaryConditions/Pressure2.cu b/src/GksGpu/BoundaryConditions/Pressure2.cu
new file mode 100644
index 0000000000000000000000000000000000000000..3aa810b3cec5ff0a5db4d12dab30b36649723b89
--- /dev/null
+++ b/src/GksGpu/BoundaryConditions/Pressure2.cu
@@ -0,0 +1,177 @@
+#include "Pressure2.h"
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include "Core/PointerDefinitions.h"
+#include "Core/RealConstants.h"
+
+#include "DataBase/DataBase.h"
+#include "DataBase/DataBaseStruct.h"
+
+#include "Definitions/MemoryAccessPattern.h"
+#include "Definitions/PassiveScalar.h"
+
+#include "FlowStateData/FlowStateData.cuh"
+#include "FlowStateData/FlowStateDataConversion.cuh"
+#include "FlowStateData/AccessDeviceData.cuh"
+
+#include "CudaUtility/CudaRunKernel.hpp"
+
+//////////////////////////////////////////////////////////////////////////
+
+__global__ void boundaryConditionKernel ( const DataBaseStruct dataBase,
+ const Pressure2Struct boundaryCondition,
+ const Parameters parameters,
+ const uint startIndex,
+ const uint numberOfEntities );
+
+__host__ __device__ inline void boundaryConditionFunction( const DataBaseStruct& dataBase,
+ const Pressure2Struct& boundaryCondition,
+ const Parameters& parameters,
+ const uint startIndex,
+ const uint index );
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+void Pressure2::runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
+ const Parameters parameters,
+ const uint level)
+{
+ CudaUtility::CudaGrid grid( this->numberOfCellsPerLevel[ level ], 32 );
+
+ runKernel( boundaryConditionKernel,
+ boundaryConditionFunction,
+ dataBase->getDeviceType(), grid,
+ dataBase->toStruct(),
+ this->toStruct(),
+ parameters,
+ this->startOfCellsPerLevel[ level ] );
+
+ cudaDeviceSynchronize();
+
+ getLastCudaError("Pressure::runBoundaryConditionKernel( const SPtr<DataBase> dataBase, const Parameters parameters, const uint level )");
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+__global__ void boundaryConditionKernel(const DataBaseStruct dataBase,
+ const Pressure2Struct boundaryCondition,
+ const Parameters parameters,
+ const uint startIndex,
+ const uint numberOfEntities)
+{
+ uint index = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if( index >= numberOfEntities ) return;
+
+ boundaryConditionFunction( dataBase, boundaryCondition, parameters, startIndex, index );
+}
+
+__host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct& dataBase,
+ const Pressure2Struct& boundaryCondition,
+ const Parameters& parameters,
+ const uint startIndex,
+ const uint index)
+{
+ uint ghostCellIdx = boundaryCondition.ghostCells [ startIndex + index ];
+ uint domainCellIdx = boundaryCondition.domainCells[ startIndex + index ];
+ uint secondCellIdx = boundaryCondition.secondCells[ startIndex + index ];
+
+ PrimitiveVariables ghostCellPrim;
+ {
+ PrimitiveVariables domainCellPrim;
+ PrimitiveVariables secondCellPrim;
+
+ {
+ ConservedVariables domainCellData;
+ readCellData( domainCellIdx, dataBase, domainCellData );
+ domainCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
+
+ ConservedVariables secondCellData;
+ if( secondCellIdx != INVALID_INDEX ){
+ readCellData( secondCellIdx, dataBase, secondCellData );
+ secondCellPrim = toPrimitiveVariables( secondCellData, parameters.K );
+ }
+ }
+
+ //ghostCellPrim.rho = two * domainCellPrim.rho - secondCellPrim.rho;
+ ghostCellPrim.U = two * domainCellPrim.U - secondCellPrim.U;
+ ghostCellPrim.V = two * domainCellPrim.V - secondCellPrim.V;
+ ghostCellPrim.W = two * domainCellPrim.W - secondCellPrim.W;
+ //ghostCellPrim.lambda = two * domainCellPrim.lambda - secondCellPrim.lambda;
+ ghostCellPrim.lambda = domainCellPrim.lambda;
+ #ifdef USE_PASSIVE_SCALAR
+ ghostCellPrim.S_1 = two * domainCellPrim.S_1 - secondCellPrim.S_1;
+ ghostCellPrim.S_2 = two * domainCellPrim.S_2 - secondCellPrim.S_2;
+ #endif // USE_PASSIVE_SCALAR
+
+
+ real rho0 = ( two * boundaryCondition.p0 * c1o2 * ( domainCellPrim.lambda + ghostCellPrim.lambda ) );
+ ghostCellPrim.rho = two * rho0 - domainCellPrim.rho;
+
+ //real lambda0 = ( c1o2 * ( domainCellPrim.rho + ghostCellPrim.rho ) * c1o2 / boundaryCondition.p0 );
+ //ghostCellPrim.lambda = two * lambda0 - domainCellPrim.lambda;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ real xGhostCell = dataBase.cellCenter[VEC_X(ghostCellIdx, dataBase.numberOfCells)];
+ real yGhostCell = dataBase.cellCenter[VEC_Y(ghostCellIdx, dataBase.numberOfCells)];
+ real zGhostCell = dataBase.cellCenter[VEC_Z(ghostCellIdx, dataBase.numberOfCells)];
+
+ real xDomainCell = dataBase.cellCenter[VEC_X(domainCellIdx, dataBase.numberOfCells)];
+ real yDomainCell = dataBase.cellCenter[VEC_Y(domainCellIdx, dataBase.numberOfCells)];
+ real zDomainCell = dataBase.cellCenter[VEC_Z(domainCellIdx, dataBase.numberOfCells)];
+
+ real dx = xGhostCell - xDomainCell;
+ real dy = yGhostCell - yDomainCell;
+ real dz = zGhostCell - zDomainCell;
+
+ real sign = domainCellPrim.U * dx
+ + domainCellPrim.V * dy
+ + domainCellPrim.W * dz;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ if( sign < zero )
+ {
+ ghostCellPrim.U = - domainCellPrim.U;
+ ghostCellPrim.V = - domainCellPrim.V;
+ ghostCellPrim.W = - domainCellPrim.W;
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+
+ {
+ ConservedVariables ghostCons = toConservedVariables( ghostCellPrim, parameters.K );
+
+ writeCellData( ghostCellIdx, dataBase, ghostCons );
+ }
+}
+
+Pressure2::Pressure2(SPtr<DataBase> dataBase, real p0)
+ : BoundaryCondition( dataBase )
+{
+ this->p0 = p0;
+}
+
+bool Pressure2::isWall()
+{
+ return false;
+}
+
+bool Pressure2::secondCellsNeeded()
+{
+ return true;
+}
+
diff --git a/src/GksGpu/BoundaryConditions/Pressure2.h b/src/GksGpu/BoundaryConditions/Pressure2.h
new file mode 100644
index 0000000000000000000000000000000000000000..59bede385b456c4d02981ba55f994e1a2c860166
--- /dev/null
+++ b/src/GksGpu/BoundaryConditions/Pressure2.h
@@ -0,0 +1,62 @@
+#ifndef Pressure2_CUH
+#define Pressure2_CUH
+
+#include <memory>
+
+#include "VirtualFluidsDefinitions.h"
+
+#include "Core/PointerDefinitions.h"
+#include "Core/DataTypes.h"
+#include "Core/VectorTypes.h"
+
+#include "BoundaryConditions/BoundaryCondition.h"
+
+//struct IsothermalWallStruct : virtual public BoundaryConditionStruct
+//{
+// Vec3 velocity;
+// real lambda;
+// real S;
+//};
+
+struct Pressure2Struct
+{
+ uint numberOfCells;
+
+ uint* ghostCells;
+ uint* domainCells;
+ uint* secondCells;
+
+ real p0;
+};
+
+struct VF_PUBLIC Pressure2 : public BoundaryCondition
+{
+ real p0;
+
+ Pressure2( SPtr<DataBase> dataBase, real p0 );
+
+ virtual bool isWall() override;
+
+ virtual bool secondCellsNeeded() override;
+
+ virtual void runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
+ const Parameters parameters,
+ const uint level) override;
+
+ Pressure2Struct toStruct()
+ {
+ Pressure2Struct boundaryCondition;
+
+ boundaryCondition.numberOfCells = this->numberOfCells;
+
+ boundaryCondition.ghostCells = this->ghostCells;
+ boundaryCondition.domainCells = this->domainCells;
+ boundaryCondition.secondCells = this->secondCells;
+
+ boundaryCondition.p0 = this->p0;
+
+ return boundaryCondition;
+ }
+};
+
+#endif
diff --git a/src/GksGpu/BoundaryConditions/Symmetry.cu b/src/GksGpu/BoundaryConditions/Symmetry.cu
new file mode 100644
index 0000000000000000000000000000000000000000..152b4239574973b7d4abaae762c99002ab3008bc
--- /dev/null
+++ b/src/GksGpu/BoundaryConditions/Symmetry.cu
@@ -0,0 +1,124 @@
+#include "Symmetry.h"
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include "Core/PointerDefinitions.h"
+#include "Core/RealConstants.h"
+
+#include "DataBase/DataBase.h"
+#include "DataBase/DataBaseStruct.h"
+
+#include "Definitions/MemoryAccessPattern.h"
+#include "Definitions/PassiveScalar.h"
+
+#include "FlowStateData/FlowStateData.cuh"
+#include "FlowStateData/FlowStateDataConversion.cuh"
+#include "FlowStateData/AccessDeviceData.cuh"
+
+#include "CudaUtility/CudaRunKernel.hpp"
+
+//////////////////////////////////////////////////////////////////////////
+
+__global__ void boundaryConditionKernel ( const DataBaseStruct dataBase,
+ const SymmetryStruct boundaryCondition,
+ const Parameters parameters,
+ const uint startIndex,
+ const uint numberOfEntities );
+
+__host__ __device__ inline void boundaryConditionFunction( const DataBaseStruct& dataBase,
+ const SymmetryStruct& boundaryCondition,
+ const Parameters& parameters,
+ const uint startIndex,
+ const uint index );
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+void Symmetry::runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
+ const Parameters parameters,
+ const uint level)
+{
+ CudaUtility::CudaGrid grid( this->numberOfCellsPerLevel[ level ], 32 );
+
+ runKernel( boundaryConditionKernel,
+ boundaryConditionFunction,
+ dataBase->getDeviceType(), grid,
+ dataBase->toStruct(),
+ this->toStruct(),
+ parameters,
+ this->startOfCellsPerLevel[ level ] );
+
+ getLastCudaError("IsothermalWall::runBoundaryConditionKernel( const SPtr<DataBase> dataBase, const Parameters parameters, const uint level )");
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+__global__ void boundaryConditionKernel(const DataBaseStruct dataBase,
+ const SymmetryStruct boundaryCondition,
+ const Parameters parameters,
+ const uint startIndex,
+ const uint numberOfEntities)
+{
+ uint index = blockIdx.x * blockDim.x + threadIdx.x;
+
+ if( index >= numberOfEntities ) return;
+
+ boundaryConditionFunction( dataBase, boundaryCondition, parameters, startIndex, index );
+}
+
+__host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct& dataBase,
+ const SymmetryStruct& boundaryCondition,
+ const Parameters& parameters,
+ const uint startIndex,
+ const uint index)
+{
+ uint ghostCellIdx = boundaryCondition.ghostCells [ startIndex + index ];
+ uint domainCellIdx = boundaryCondition.domainCells[ startIndex + index ];
+
+ //////////////////////////////////////////////////////////////////////////
+
+ ConservedVariables domainCellData;
+ readCellData ( domainCellIdx, dataBase, domainCellData );
+
+ //////////////////////////////////////////////////////////////////////////
+
+ PrimitiveVariables domainCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
+ PrimitiveVariables ghostCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
+
+ //////////////////////////////////////////////////////////////////////////
+
+ if( boundaryCondition.direction == 'x' ) ghostCellPrim.U = - domainCellPrim.U;
+ if( boundaryCondition.direction == 'y' ) ghostCellPrim.V = - domainCellPrim.V;
+ if( boundaryCondition.direction == 'z' ) ghostCellPrim.W = - domainCellPrim.W;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ ConservedVariables ghostCellData = toConservedVariables(ghostCellPrim, parameters.K);
+
+ writeCellData( ghostCellIdx , dataBase, ghostCellData );
+}
+
+Symmetry::Symmetry(SPtr<DataBase> dataBase, char direction)
+ : BoundaryCondition( dataBase )
+{
+ this->direction = direction;
+}
+
+bool Symmetry::isWall()
+{
+ return true;
+}
+
+bool Symmetry::secondCellsNeeded()
+{
+ return false;
+}
+
diff --git a/src/GksGpu/BoundaryConditions/Symmetry.h b/src/GksGpu/BoundaryConditions/Symmetry.h
new file mode 100644
index 0000000000000000000000000000000000000000..165ad040a2967d9c2184419c32c10fa25be1cba2
--- /dev/null
+++ b/src/GksGpu/BoundaryConditions/Symmetry.h
@@ -0,0 +1,55 @@
+#ifndef Symmetry_CUH
+#define Symmetry_CUH
+
+#include <memory>
+
+#include "VirtualFluidsDefinitions.h"
+
+#include "Core/PointerDefinitions.h"
+#include "Core/DataTypes.h"
+#include "Core/VectorTypes.h"
+
+#include "BoundaryConditions/BoundaryCondition.h"
+
+struct SymmetryStruct
+{
+ uint numberOfCells;
+
+ uint* ghostCells;
+ uint* domainCells;
+ uint* secondCells;
+
+ char direction;
+};
+
+struct VF_PUBLIC Symmetry : public BoundaryCondition //, public IsothermalWallStruct
+{
+ char direction;
+
+ Symmetry( SPtr<DataBase> dataBase, char direction );
+
+ virtual bool isWall() override;
+
+ virtual bool secondCellsNeeded() override;
+
+ virtual void runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
+ const Parameters parameters,
+ const uint level) override;
+
+ SymmetryStruct toStruct()
+ {
+ SymmetryStruct boundaryCondition;
+
+ boundaryCondition.numberOfCells = this->numberOfCells;
+
+ boundaryCondition.ghostCells = this->ghostCells;
+ boundaryCondition.domainCells = this->domainCells;
+ boundaryCondition.secondCells = this->secondCells;
+
+ boundaryCondition.direction = this->direction;
+
+ return boundaryCondition;
+ }
+};
+
+#endif
diff --git a/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp b/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
index 82d584cc9650e322d86ee7539f9aba0d15d94018..051d9515219117cdbef94bc27c866bc68d2bf7f7 100644
--- a/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
+++ b/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
@@ -131,6 +131,7 @@ void DataBaseAllocatorCPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
dataBase->faceCenter[ VEC_X( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.x;
dataBase->faceCenter[ VEC_Y( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.y;
+ dataBase->faceCenter[ VEC_Z( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.z;
dataBase->faceOrientation[ faceIdx ] = adapter.faces[ faceIdx ].orientation;
}
diff --git a/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp b/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
index 63012398d3d2676430220d737f5130feefff1876..d8e33597bfad01ff293a10a8b858b5dd63a79b5b 100644
--- a/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
+++ b/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
@@ -152,6 +152,7 @@ void DataBaseAllocatorGPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
faceCenterBuffer[ VEC_X( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.x;
faceCenterBuffer[ VEC_Y( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.y;
+ faceCenterBuffer[ VEC_Z( faceIdx, dataBase->numberOfFaces ) ] = adapter.faces[ faceIdx ].faceCenter.z;
faceOrientationBuffer[ faceIdx ] = adapter.faces[ faceIdx ].orientation;
}
diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index 2a596b0df1ef71d8ad6855f3e621e069b82b2d13..3b72682fe6ed5cfb78d7e450d04c841c6192478b 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -122,13 +122,13 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
real muNew = parameters.mu;
- if( fabsf(x) > c3o2 )
+ if( fabsf(x) > three )
{
- muNew += ( fabsf(x) - c3o2 ) * ten * parameters.mu;
+ muNew += ( fabsf(x) - three ) * ten * parameters.mu;
}
- if( fabsf(z) > three )
+ if( fabsf(z) > seven )
{
- muNew += ( fabs(z) - three ) * ten * parameters.mu;
+ muNew += ( fabs(z) - seven ) * ten * parameters.mu;
}
//parameters.Pr = muNew / parameters.mu;
diff --git a/targets/apps/GKS/PoolFire/PoolFire.cpp b/targets/apps/GKS/PoolFire/PoolFire.cpp
index 7cf60ea326c60cd8ecb193873846e3072b8e8028..10c000c6bd7dff1987707d1798a12ed954459abf 100644
--- a/targets/apps/GKS/PoolFire/PoolFire.cpp
+++ b/targets/apps/GKS/PoolFire/PoolFire.cpp
@@ -37,11 +37,13 @@
#include "GksGpu/BoundaryConditions/BoundaryCondition.h"
#include "GksGpu/BoundaryConditions/IsothermalWall.h"
#include "GksGpu/BoundaryConditions/Periodic.h"
-#include "GksGpu/BoundaryConditions/Pressure.h"
+#include "GksGpu/BoundaryConditions/Pressure2.h"
#include "GksGpu/BoundaryConditions/AdiabaticWall.h"
#include "GksGpu/BoundaryConditions/PassiveScalarDiriclet.h"
#include "GksGpu/BoundaryConditions/InflowComplete.h"
#include "GksGpu/BoundaryConditions/Open.h"
+#include "GksGpu/BoundaryConditions/Extrapolation.h"
+#include "GksGpu/BoundaryConditions/Symmetry.h"
#include "GksGpu/Interface/Interface.h"
#include "GksGpu/TimeStepping/NestedTimeStep.h"
@@ -50,9 +52,11 @@
#include "GksGpu/Analyzer/ConvergenceAnalyzer.h"
#include "GksGpu/Analyzer/TurbulenceAnalyzer.h"
+#include "GksGpu/Restart/Restart.h"
+
#include "GksGpu/CudaUtility/CudaUtility.h"
-void thermalCavity( std::string path, std::string simulationName )
+void thermalCavity( std::string path, std::string simulationName, uint restartIter )
{
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -128,7 +132,7 @@ void thermalCavity( std::string path, std::string simulationName )
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- bool threeDimensional = false;
+ bool threeDimensional = true;
if( threeDimensional )
{
@@ -191,7 +195,8 @@ void thermalCavity( std::string path, std::string simulationName )
//meshAdapter.writeMeshFaceVTK( path + "grid/MeshFaces.vtk" );
- meshAdapter.findPeriodicBoundaryNeighbors();
+ if( !threeDimensional )
+ meshAdapter.findPeriodicBoundaryNeighbors();
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -206,14 +211,23 @@ void thermalCavity( std::string path, std::string simulationName )
SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase, prim );
SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase, prim );
- SPtr<BoundaryCondition> bcMX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
- SPtr<BoundaryCondition> bcPX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
+ //SPtr<BoundaryCondition> bcMX = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
+ //SPtr<BoundaryCondition> bcPX = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
+ //SPtr<BoundaryCondition> bcMX = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
+ //SPtr<BoundaryCondition> bcPX = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
bcMX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L; } );
bcPX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L; } );
- bcMX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L && center.z > 0.5; } );
- bcPX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L && center.z > 0.5; } );
+ //SPtr<BoundaryCondition> bcMX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
+ //SPtr<BoundaryCondition> bcPX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
+ SPtr<BoundaryCondition> bcMX_2 = std::make_shared<Symmetry>( dataBase, 'x' );
+ SPtr<BoundaryCondition> bcPX_2 = std::make_shared<Symmetry>( dataBase, 'x' );
+ //SPtr<BoundaryCondition> bcMX_2 = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
+ //SPtr<BoundaryCondition> bcPX_2 = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
+
+ bcMX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L && center.z > 1.0; } );
+ bcPX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L && center.z > 1.0; } );
//////////////////////////////////////////////////////////////////////////
@@ -222,8 +236,10 @@ void thermalCavity( std::string path, std::string simulationName )
if( threeDimensional )
{
- bcMY = std::make_shared<Open>( dataBase, prim );
- bcPY = std::make_shared<Open>( dataBase, prim );
+ //bcMY = std::make_shared<Open>( dataBase, prim );
+ //bcPY = std::make_shared<Open>( dataBase, prim );
+ bcMY = std::make_shared<Symmetry>( dataBase, 'y' );
+ bcPY = std::make_shared<Symmetry>( dataBase, 'y' );
bcMY->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.y < -0.5*L; } );
bcPY->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.y > 0.5*L; } );
@@ -244,7 +260,9 @@ void thermalCavity( std::string path, std::string simulationName )
//SPtr<BoundaryCondition> bcMZ = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, 0.0, prim.lambda, 0.0, 0.0) );
//SPtr<BoundaryCondition> bcMZ = std::make_shared<Open>( dataBase );
- SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
+ //SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
+ //SPtr<BoundaryCondition> bcPZ = std::make_shared<Extrapolation>( dataBase );
+ SPtr<BoundaryCondition> bcPZ = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0; } );
bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z > H ; } );
@@ -266,8 +284,8 @@ void thermalCavity( std::string path, std::string simulationName )
dataBase->boundaryConditions.push_back( bcMX );
dataBase->boundaryConditions.push_back( bcPX );
- //dataBase->boundaryConditions.push_back( bcMX_2 );
- //dataBase->boundaryConditions.push_back( bcPX_2 );
+ dataBase->boundaryConditions.push_back( bcMX_2 );
+ dataBase->boundaryConditions.push_back( bcPX_2 );
dataBase->boundaryConditions.push_back( bcMY );
dataBase->boundaryConditions.push_back( bcPY );
@@ -282,24 +300,37 @@ void thermalCavity( std::string path, std::string simulationName )
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ uint startIter = 0;
+
dataBase->setMesh( meshAdapter );
CudaUtility::printCudaMemoryUsage();
+
+ if( restartIter == INVALID_INDEX )
+ {
+ Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> ConservedVariables {
- Initializer::interpret(dataBase, [&] ( Vec3 cellCenter ) -> ConservedVariables{
+ PrimitiveVariables primLocal = prim;
- PrimitiveVariables primLocal = prim;
-
- //primLocal.rho = rho * std::exp( - ( 2.0 * g * H * prim.lambda ) * cellCenter.z / H );
+ //primLocal.rho = rho * std::exp( - ( 2.0 * g * H * prim.lambda ) * cellCenter.z / H );
- real r = sqrt( cellCenter.x * cellCenter.x + cellCenter.y * cellCenter.y /*+ cellCenter.z * cellCenter.z*/ );
+ real r = sqrt(cellCenter.x * cellCenter.x + cellCenter.y * cellCenter.y /*+ cellCenter.z * cellCenter.z*/);
- //if( r < 0.6 ) primLocal.S_1 = 1.0;
+ //if( r < 0.6 ) primLocal.S_1 = 1.0;
- //if( r < 0.5 ) prim.lambda /= (two - four*r*r);
+ //if( r < 0.5 ) prim.lambda /= (two - four*r*r);
- return toConservedVariables( primLocal, parameters.K );
- });
+ return toConservedVariables(primLocal, parameters.K);
+ });
+
+ writeVtkXML( dataBase, parameters, 0, path + simulationName + "_0" );
+ }
+ else
+ {
+ Restart::readRestart( dataBase, path + simulationName + "_" + std::to_string( restartIter ), startIter );
+
+ writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( restartIter ) + "_restart" );
+ }
dataBase->copyDataHostToDevice();
@@ -311,8 +342,6 @@ void thermalCavity( std::string path, std::string simulationName )
dataBase->copyDataDeviceToHost();
- writeVtkXML( dataBase, parameters, 0, path + simulationName + "_0" );
-
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -328,7 +357,7 @@ void thermalCavity( std::string path, std::string simulationName )
cupsAnalyzer.start();
- for( uint iter = 1; iter <= 100000000; iter++ )
+ for( uint iter = startIter + 1; iter <= 100000000; iter++ )
{
if( iter < 20000 )
{
@@ -368,6 +397,11 @@ void thermalCavity( std::string path, std::string simulationName )
writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
}
+ if( iter % 10000 == 0 )
+ {
+ Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ), iter );
+ }
+
//turbulenceAnalyzer->run( iter, parameters );
}
@@ -404,7 +438,12 @@ int main( int argc, char* argv[])
try
{
- thermalCavity( path, simulationName );
+ uint restartIter = INVALID_INDEX;
+ //uint restartIter = 60000;
+
+ if( argc > 1 ) restartIter = atoi( argv[1] );
+
+ thermalCavity( path, simulationName, restartIter );
}
catch (const std::exception& e)
{