seemingly plausible state tested in 2D and 3D, crashes due to insufficient...

seemingly plausible state tested in 2D and 3D, crashes due to insufficient lateral boundary conditions

src/GksGpu/BoundaryConditions/InflowComplete.cu

@@ -99,7 +99,7 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
     //    ghostCellPrim.U      = two * boundaryCondition.prim.U      - domainCellPrim.U;
     //    ghostCellPrim.V      = two * boundaryCondition.prim.V      - domainCellPrim.V;
     //    ghostCellPrim.W      = two * boundaryCondition.prim.W      - domainCellPrim.W;
-    //    ghostCellPrim.lambda = /*two * boundaryCondition.prim.lambda -*/ domainCellPrim.lambda;
+        ghostCellPrim.lambda = /*two * boundaryCondition.prim.lambda -*/ domainCellPrim.lambda;
     //#ifdef USE_PASSIVE_SCALAR
     //    ghostCellPrim.S_1    = two * boundaryCondition.prim.S_1    - domainCellPrim.S_1;
     //    ghostCellPrim.S_2    = two * boundaryCondition.prim.S_2    - domainCellPrim.S_2;
@@ -109,7 +109,7 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
         ghostCellPrim.U      = two * boundaryCondition.prim.U - domainCellPrim.U;
         ghostCellPrim.V      = two * boundaryCondition.prim.V - domainCellPrim.V;
         ghostCellPrim.W      = boundaryCondition.prim.W;
-        ghostCellPrim.lambda = boundaryCondition.prim.lambda;
+        //ghostCellPrim.lambda = boundaryCondition.prim.lambda;
     #ifdef USE_PASSIVE_SCALAR
         ghostCellPrim.S_1    = boundaryCondition.prim.S_1;
         ghostCellPrim.S_2    = boundaryCondition.prim.S_2;
@@ -121,7 +121,7 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
         real r = sqrt( y*y + x*x );
 
         ghostCellPrim.W   *= (one - four*r*r);
-        ghostCellPrim.S_1 *=       (one - four*r*r);
+        ghostCellPrim.S_1 *= (one - four*r*r);
         ghostCellPrim.S_2 = one - ghostCellPrim.S_1;
     }
 

src/GksGpu/BoundaryConditions/Open.cu

@@ -17,6 +17,7 @@
 #include "Definitions/PassiveScalar.h"
 
 #include "FlowStateData/FlowStateData.cuh"
+#include "FlowStateData/FlowStateDataConversion.cuh"
 #include "FlowStateData/AccessDeviceData.cuh"
 
 #include "CudaUtility/CudaRunKernel.hpp"
@@ -82,14 +83,47 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
     uint ghostCellIdx  = boundaryCondition.ghostCells [ startIndex + index ];
     uint domainCellIdx = boundaryCondition.domainCells[ startIndex + index ];
 
-    ConservedVariables domainCellData;
+    ConservedVariables domainCellData, ghostCellData;
     readCellData ( domainCellIdx, dataBase, domainCellData );
-    writeCellData( ghostCellIdx , dataBase, domainCellData );
+
+    PrimitiveVariables domainCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
+    
+    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 )
+        ghostCellData = domainCellData;
+    else
+    {
+        PrimitiveVariables ghostCellPrim  = boundaryCondition.prim;
+
+        ghostCellPrim.U = domainCellPrim.U;
+        ghostCellPrim.V = domainCellPrim.V;
+        ghostCellPrim.W = domainCellPrim.W;
+
+        ghostCellData = toConservedVariables(ghostCellPrim, parameters.K);
+    }
+
+    writeCellData(ghostCellIdx, dataBase, ghostCellData);
 }
 
-Open::Open(SPtr<DataBase> dataBase)
+Open::Open(SPtr<DataBase> dataBase, PrimitiveVariables prim)
     : BoundaryCondition( dataBase )
 {
+    this->prim = prim;
 }
 
 bool Open::isWall()

src/GksGpu/BoundaryConditions/Open.h

@@ -9,6 +9,8 @@
 #include "Core/DataTypes.h"
 #include "Core/VectorTypes.h"
 
+#include "FlowStateData/FlowStateData.cuh"
+
 #include "BoundaryConditions/BoundaryCondition.h"
 
 struct OpenStruct
@@ -18,11 +20,15 @@ struct OpenStruct
     uint* ghostCells;
     uint* domainCells;
     uint* secondCells;
+
+    PrimitiveVariables prim;
 };
 
 struct VF_PUBLIC Open : public BoundaryCondition //, public IsothermalWallStruct
 {
-    Open( SPtr<DataBase> dataBase );
+    PrimitiveVariables prim;
+
+    Open( SPtr<DataBase> dataBase, PrimitiveVariables prim );
 
     virtual bool isWall() override;
 
@@ -42,6 +48,8 @@ struct VF_PUBLIC Open : public BoundaryCondition //, public IsothermalWallStruct
         boundaryCondition.domainCells   = this->domainCells;
         boundaryCondition.secondCells   = this->secondCells;
 
+        boundaryCondition.prim          = prim;
+
         return boundaryCondition;
     }
 };

targets/apps/GKS/PoolFire/PoolFire.cpp

@@ -128,7 +128,7 @@ void thermalCavity( std::string path, std::string simulationName )
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    bool threeDimensional = false;
+    bool threeDimensional = true;
 
     if( threeDimensional )
     {
@@ -197,8 +197,10 @@ void thermalCavity( std::string path, std::string simulationName )
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     
-    SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase );
-    SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase );
+    SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase, prim );
+    SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase, prim );
+    //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 );
 
     bcMX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L; } );
     bcPX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x >  0.5*L; } );
@@ -210,8 +212,8 @@ void thermalCavity( std::string path, std::string simulationName )
 
     if( threeDimensional )
     {
-        bcMY = std::make_shared<Open>( dataBase );
-        bcPY = std::make_shared<Open>( dataBase );
+        bcMY = std::make_shared<Open>( dataBase, prim );
+        bcPY = std::make_shared<Open>( dataBase, prim );
 
         bcMY->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.y < -0.5*L; } );
         bcPY->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.y >  0.5*L; } );
@@ -232,7 +234,7 @@ 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 );
+    SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
     
     bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0; } );
     bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z > H  ; } );
@@ -331,7 +333,7 @@ void thermalCavity( std::string path, std::string simulationName )
         TimeStepping::nestedTimeStep(dataBase, parameters, 0);
 
         if( 
-            //( iter >= 10200 && iter % 10 == 0 ) || 
+            //( iter >= 5000 && iter % 10 == 0 ) || 
             ( iter % 1000 == 0 )
           )
         {