diff --git a/src/GksGpu/BoundaryConditions/BoundaryCondition.cpp b/src/GksGpu/BoundaryConditions/BoundaryCondition.cpp
index e3bb858153f14f5b61bbca3a03de0cfb94469f40..0fc85a98d243ca6ff27aea1be052b9df3f6a4d27 100644
--- a/src/GksGpu/BoundaryConditions/BoundaryCondition.cpp
+++ b/src/GksGpu/BoundaryConditions/BoundaryCondition.cpp
@@ -54,7 +54,13 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
 
             if( cell.type != STOPPER_OUT_OF_GRID && cell.type != STOPPER_OUT_OF_GRID_BOUNDARY ) continue;
 
-            for( uint idx = 0; idx < 27; idx++ )
+            // look in all directions
+            uint maximalSearchDirection = 27;
+
+            // in case of Flux BC look only at face neighbors
+            if( this->isFluxBC() ) maximalSearchDirection = 6;
+
+            for( uint idx = 0; idx < maximalSearchDirection; idx++ )
             {
                 uint neighborCellIdx = cell.cellToCell[ idx ];
 
@@ -65,7 +71,7 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
                 bool neighborCellIsFluid = neighborCell.type != STOPPER_OUT_OF_GRID && 
                                            neighborCell.type != STOPPER_OUT_OF_GRID_BOUNDARY;
 
-                bool neighborCellIsValidGhostCell = allowGhostCells && !boundaryFinder( neighborCell.cellCenter );
+                bool neighborCellIsValidGhostCell = !this->isFluxBC() && allowGhostCells && !boundaryFinder( neighborCell.cellCenter );
 
                 if( neighborCellIsFluid || neighborCellIsValidGhostCell )
                 {
@@ -79,7 +85,8 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
                         secondCells.push_back( neighborCell.cellToCell[ idx ] );
                     }
 
-                    cell.isWall = this->isWall();
+                    cell.isWall   = this->isWall();
+                    cell.isFluxBC = this->isFluxBC();
 
                     break;
                 }
@@ -100,6 +107,11 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
     this->myAllocator->allocateMemory( shared_from_this(), ghostCells, domainCells, secondCells );
 }
 
+bool BoundaryCondition::isFluxBC()
+{
+    return false;
+}
+
 bool BoundaryCondition::secondCellsNeeded()
 {
     return false;
diff --git a/src/GksGpu/BoundaryConditions/BoundaryCondition.h b/src/GksGpu/BoundaryConditions/BoundaryCondition.h
index bb4e33a37ee9bf28510cb4de15f51cba28583cca..6e53558c332dba0de56debb1cd86d54850aca27d 100644
--- a/src/GksGpu/BoundaryConditions/BoundaryCondition.h
+++ b/src/GksGpu/BoundaryConditions/BoundaryCondition.h
@@ -43,6 +43,8 @@ struct VF_PUBLIC BoundaryCondition : virtual public BoundaryConditionStruct, pub
 
     virtual bool isWall() = 0;
 
+    virtual bool isFluxBC();
+
     virtual bool secondCellsNeeded();
 
     virtual void runBoundaryConditionKernel( const SPtr<DataBase> dataBase,
diff --git a/src/GksGpu/BoundaryConditions/InflowComplete.cu b/src/GksGpu/BoundaryConditions/InflowComplete.cu
index b12fab5556fe8caf8a22df1d7030d5316c7fffa1..543e0046fcffb5969d6c356cb0b5e1415e908d0f 100644
--- a/src/GksGpu/BoundaryConditions/InflowComplete.cu
+++ b/src/GksGpu/BoundaryConditions/InflowComplete.cu
@@ -2,6 +2,7 @@
 
 #define _USE_MATH_DEFINES
 #include <math.h>
+#include <iostream>
 
 #include <cuda.h>
 #include <cuda_runtime.h>
@@ -20,6 +21,12 @@
 #include "FlowStateData/FlowStateDataConversion.cuh"
 #include "FlowStateData/AccessDeviceData.cuh"
 
+#include "FluxComputation/Moments.cuh"
+#include "FluxComputation/ApplyFlux.cuh"
+#include "FluxComputation/Transformation.cuh"
+#include "FluxComputation/AssembleFlux.cuh"
+#include "FluxComputation/ExpansionCoefficients.cuh"
+
 #include "CudaUtility/CudaRunKernel.hpp"
 
 //////////////////////////////////////////////////////////////////////////
@@ -85,53 +92,207 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
     uint ghostCellIdx  = boundaryCondition.ghostCells [ startIndex + index ];
     uint domainCellIdx = boundaryCondition.domainCells[ startIndex + index ];
 
-    PrimitiveVariables ghostCellPrim;
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    if( true )
     {
-        PrimitiveVariables domainCellPrim;
+        PrimitiveVariables ghostCellPrim;
+        {
+            PrimitiveVariables domainCellPrim;
+
+            {
+                ConservedVariables domainCellData;
+                readCellData(domainCellIdx, dataBase, domainCellData);
+                domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
+            }
+
+            //    ghostCellPrim.rho    = two * boundaryCondition.prim.rho    - domainCellPrim.rho;
+            //    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;
+            //#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;
+            //#endif // USE_PASSIVE_SCALAR
+
+            ghostCellPrim.rho = boundaryCondition.prim.rho;
+            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.W      = boundaryCondition.prim.W;
+            //ghostCellPrim.lambda = boundaryCondition.prim.lambda;
+#ifdef USE_PASSIVE_SCALAR
+            ghostCellPrim.S_1 = boundaryCondition.prim.S_1;
+            ghostCellPrim.S_2 = boundaryCondition.prim.S_2;
+#endif // USE_PASSIVE_SCALAR
+
+            real y = dataBase.cellCenter[VEC_Y(ghostCellIdx, dataBase.numberOfCells)];
+            real x = dataBase.cellCenter[VEC_X(ghostCellIdx, dataBase.numberOfCells)];
+
+            real r = sqrt(y*y + x*x);
+
+            ghostCellPrim.W *= (one - four*r*r);
+#ifdef USE_PASSIVE_SCALAR
+            ghostCellPrim.S_1 *= (one - four*r*r);
+            ghostCellPrim.S_2 = boundaryCondition.prim.S_2 - ghostCellPrim.S_1;
+#endif // USE_PASSIVE_SCALAR
+        }
+
+        {
+            ConservedVariables ghostCons = toConservedVariables(ghostCellPrim, parameters.K);
+
+            writeCellData(ghostCellIdx, dataBase, ghostCons);
+        }
+    }
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+    if( false )
+    {
+        PrimitiveVariables domainCellPrim;
         {
             ConservedVariables domainCellData;
-            readCellData( domainCellIdx, dataBase, domainCellData );
-            domainCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
+            readCellData(domainCellIdx, dataBase, domainCellData);
+            domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
         }
 
-    //    ghostCellPrim.rho    = two * boundaryCondition.prim.rho    - domainCellPrim.rho;
-    //    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;
-    //#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;
-    //#endif // USE_PASSIVE_SCALAR
-
-        ghostCellPrim.rho    = boundaryCondition.prim.rho;
-        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;
-    #ifdef USE_PASSIVE_SCALAR
-        ghostCellPrim.S_1    = boundaryCondition.prim.S_1;
-        ghostCellPrim.S_2    = boundaryCondition.prim.S_2;
-    #endif // USE_PASSIVE_SCALAR
-
-        real y = dataBase.cellCenter[ VEC_Y(ghostCellIdx, dataBase.numberOfCells) ];
-        real x = dataBase.cellCenter[ VEC_X(ghostCellIdx, dataBase.numberOfCells) ];
-
-        real r = sqrt( y*y + x*x );
-
-        ghostCellPrim.W   *= (one - four*r*r);
-    #ifdef USE_PASSIVE_SCALAR
-        ghostCellPrim.S_1 *= (one - four*r*r);
-        ghostCellPrim.S_2 = one - ghostCellPrim.S_1;
-    #endif // USE_PASSIVE_SCALAR
+        real momentU [NUMBER_OF_MOMENTS];
+        real momentV [NUMBER_OF_MOMENTS];
+        real momentW [NUMBER_OF_MOMENTS];
+        real momentXi[NUMBER_OF_MOMENTS];
+
+        PrimitiveVariables facePrim = boundaryCondition.prim;
+
+        //facePrim.lambda = domainCellPrim.lambda;
+
+        transformGlobalToLocal( facePrim, 'z' );
+
+        computeMoments(facePrim, parameters.K, momentU, momentV, momentW, momentXi);
+
+        ConservedVariables flux;
+
+        flux.rho  = momentU[0 + 1];
+        //flux.rhoU = momentU[1 + 1];
+
+        //flux.rhoE = c1o2 * ( momentU[2 + 1]
+        //                   + momentU[0 + 1] * momentV [2]
+        //                   + momentU[0 + 1] * momentW [2]
+        //                   + momentU[0 + 1] * momentXi[2] );
+
+        flux.rhoE = momentU[0 + 1] * c1o4 * ( parameters.K + five ) / boundaryCondition.prim.lambda;
+
+        //////////////////////////////////////////////////////////////////////////
+
+#ifdef USE_PASSIVE_SCALAR
+        flux.rhoS_1 = flux.rho * boundaryCondition.prim.S_1;
+        flux.rhoS_2 = flux.rho * boundaryCondition.prim.S_2;
+#endif // USE_PASSIVE_SCALAR
+
+        flux   = ( parameters.dt * parameters.dx * parameters.dx * facePrim.rho ) * flux;
+
+        transformLocalToGlobal( flux, 'z' );
+
+        applyFluxToPosCell(dataBase, domainCellIdx, flux, 'z', parameters.dt);
+
+        return;
     }
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    if( false )
     {
-        ConservedVariables ghostCons = toConservedVariables( ghostCellPrim, parameters.K );
+        PrimitiveVariables domainCellPrim;
+        {
+            ConservedVariables domainCellData;
+            readCellData(domainCellIdx, dataBase, domainCellData);
+            domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
+        }    
+
+        PrimitiveVariables facePrim = boundaryCondition.prim;
+
+        //////////////////////////////////////////////////////////////////////////
+
+        real ax[LENGTH_CELL_DATA];
+        real ay[LENGTH_CELL_DATA];
+        real az[LENGTH_CELL_DATA];
+        real at[LENGTH_CELL_DATA];
+
+    #pragma unroll
+        for( uint i = 0; i < LENGTH_CELL_DATA; i++ )
+        { 
+            ax[i] = zero; 
+            ay[i] = zero; 
+            az[i] = zero; 
+            at[i] = zero;
+        }
+        
+        {
+            ConservedVariables gradN, gradT1, gradT2;
 
-        writeCellData( ghostCellIdx, dataBase, ghostCons );
+            transformGlobalToLocal( gradN , 'z' );
+            transformGlobalToLocal( gradT1, 'z' );
+            transformGlobalToLocal( gradT2, 'z' );
+
+            transformGlobalToLocal( facePrim, 'z' );
+
+            computeExpansionCoefficients(facePrim, gradN , parameters.K, ax);
+            computeExpansionCoefficients(facePrim, gradT1, parameters.K, ay);
+            computeExpansionCoefficients(facePrim, gradT2, parameters.K, az);
+        }
+
+        //////////////////////////////////////////////////////////////////////////
+
+        {
+            ConservedVariables flux;
+            {
+                real momentU [ NUMBER_OF_MOMENTS ]; 
+                real momentV [ NUMBER_OF_MOMENTS ]; 
+                real momentW [ NUMBER_OF_MOMENTS ]; 
+                real momentXi[ NUMBER_OF_MOMENTS ];
+
+                computeMoments( facePrim, parameters.K, momentU, momentV, momentW, momentXi );
+
+                Vec3 force = parameters.force;
+
+                transformGlobalToLocal(force, 'z');
+
+                {
+                    ConservedVariables timeGrad;
+                    computeTimeDerivative( facePrim, 
+                                           momentU, 
+                                           momentV, 
+                                           momentW, 
+                                           momentXi, 
+                                           ax, ay, az,
+                                           force,
+                                           timeGrad );
+
+                    computeExpansionCoefficients( facePrim, timeGrad, parameters.K, at );
+                }
+                {
+                    real timeCoefficients[4];
+                    computeTimeCoefficients( facePrim, parameters, timeCoefficients );
+
+                    real heatFlux;
+                    assembleFlux( facePrim, 
+                                  momentU, momentV, momentW, momentXi,
+                                  ax, ay, az, at, 
+                                  timeCoefficients, 
+                                  parameters,
+                                  force,
+                                  flux,
+                                  heatFlux );
+
+                    transformLocalToGlobal( flux, 'z' );
+                }
+            }
+
+            applyFluxToPosCell(dataBase, domainCellIdx, flux, 'z', parameters.dt);
+            applyFluxToNegCell(dataBase, ghostCellIdx , flux, 'z', parameters.dt);
+        }
     }
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 }
 
 InflowComplete::InflowComplete(SPtr<DataBase> dataBase, PrimitiveVariables prim)
@@ -145,6 +306,11 @@ bool InflowComplete::isWall()
     return false;
 }
 
+bool InflowComplete::isFluxBC()
+{
+    return false;
+}
+
 bool InflowComplete::secondCellsNeeded()
 {
     return false;
diff --git a/src/GksGpu/BoundaryConditions/InflowComplete.h b/src/GksGpu/BoundaryConditions/InflowComplete.h
index 1952ddabefd1b4c525b93c801893c0d6cf019863..e444bd1ba066f0915e6a4d104b4d4d5d3e673b6f 100644
--- a/src/GksGpu/BoundaryConditions/InflowComplete.h
+++ b/src/GksGpu/BoundaryConditions/InflowComplete.h
@@ -32,6 +32,8 @@ struct VF_PUBLIC InflowComplete : public BoundaryCondition //, public Isothermal
 
     virtual bool isWall() override;
 
+    virtual bool isFluxBC() override;
+
     virtual bool secondCellsNeeded() override;
 
     virtual void runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
diff --git a/src/GksGpu/BoundaryConditions/Open.cu b/src/GksGpu/BoundaryConditions/Open.cu
index 62d1b8d43e99d6e4d722cafad7166be1f502efa9..87d14e17e685aa87872c244e0164aeea8aa3303a 100644
--- a/src/GksGpu/BoundaryConditions/Open.cu
+++ b/src/GksGpu/BoundaryConditions/Open.cu
@@ -116,8 +116,8 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
 
     //////////////////////////////////////////////////////////////////////////
 
-    //if( sign > zero )
-    if( p2 > p1 )
+    if( sign < zero )
+    //if( p2 > p1 )
         ghostCellData = domainCellData;
     else
     {
@@ -139,6 +139,14 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
         ghostCellData = toConservedVariables(ghostCellPrim, parameters.K);
     }
 
+    //////////////////////////////////////////////////////////////////////////
+
+    //ghostCellData = two * domainCellData + ( - one ) * secondCellData;
+
+    //ghostCellData = domainCellData;
+
+    //////////////////////////////////////////////////////////////////////////
+
     writeCellData(ghostCellIdx, dataBase, ghostCellData);
 }
 
diff --git a/src/GksGpu/CellProperties/CellProperties.cuh b/src/GksGpu/CellProperties/CellProperties.cuh
index f50aa07edd85f89de557f8a0a566eed870024d24..83be81c19e9489148bf5143912cfaf185a79e762 100644
--- a/src/GksGpu/CellProperties/CellProperties.cuh
+++ b/src/GksGpu/CellProperties/CellProperties.cuh
@@ -14,7 +14,7 @@
 #define CELL_PROPERTIES_GHOST      (1u)
 #define CELL_PROPERTIES_WALL       (2u)
 #define CELL_PROPERTIES_FINE_GHOST (4u)
-#define CELL_PROPERTIES_3          (8u)
+#define CELL_PROPERTIES_IS_FLUX_BC (8u)
 #define CELL_PROPERTIES_4          (16u)
 #define CELL_PROPERTIES_5          (32u)
 #define CELL_PROPERTIES_6          (64u)
diff --git a/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp b/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
index 331e77d71305d0ac52d2b819f83935f1c7f11e7e..c909ae7064faf0503baee5c455f093c040da8363 100644
--- a/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
+++ b/src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
@@ -105,6 +105,8 @@ void DataBaseAllocatorCPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
         dataBase->cellProperties[ cellIdx ] = CELL_PROPERTIES_DEFAULT;
         if( adapter.cells[ cellIdx ].isWall )
             setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_WALL ); 
+        if( adapter.cells[ cellIdx ].isFluxBC )
+            setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_IS_FLUX_BC ); 
         if( adapter.cells[ cellIdx ].isGhostCell )
             setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_GHOST );
         if( adapter.cells[ cellIdx ].isFineGhostCell() )
diff --git a/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp b/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
index 4f4e022ce4d0203a053a9335b7408b24b7d9a105..27601e1d535e024066b27913406320e6a896f0fa 100644
--- a/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
+++ b/src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
@@ -126,6 +126,8 @@ void DataBaseAllocatorGPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
         cellPropertiesBuffer[ cellIdx ] = CELL_PROPERTIES_DEFAULT;
         if( adapter.cells[ cellIdx ].isWall )
             setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_WALL ); 
+        if( adapter.cells[ cellIdx ].isFluxBC )
+            setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_IS_FLUX_BC ); 
         if( adapter.cells[ cellIdx ].isGhostCell )
             setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_GHOST ); 
         if( adapter.cells[ cellIdx ].isFineGhostCell() )
diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index 58b20a8946d2fdf1ecdde05cfde5d881dd9644e0..c5522b1f1cd7b6a5691d048808c356f8c6d0d7b4 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -190,8 +190,6 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
     {
         ConservedVariables flux;
 
-        ConservedVariables faceCons;
-
         {
             real momentU [ NUMBER_OF_MOMENTS ]; 
             real momentV [ NUMBER_OF_MOMENTS ]; 
@@ -244,6 +242,10 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
             CellProperties negCellProperties = dataBase.cellProperties[ negCellIdx ];
             CellProperties posCellProperties = dataBase.cellProperties[ posCellIdx ];
 
+            //if( isCellProperties( negCellProperties, CELL_PROPERTIES_IS_FLUX_BC ) || 
+            //    isCellProperties( posCellProperties, CELL_PROPERTIES_IS_FLUX_BC ) )
+            //    return;
+
             if( isCellProperties( negCellProperties, CELL_PROPERTIES_WALL ) || 
                 isCellProperties( posCellProperties, CELL_PROPERTIES_WALL ) )
             {
diff --git a/src/GksMeshAdapter/MeshCell.cpp b/src/GksMeshAdapter/MeshCell.cpp
index ded3f193f21c88db59722487b6ef6eefce600259..04e116ea305a7e3c42a9a3085425aba848ee5034 100644
--- a/src/GksMeshAdapter/MeshCell.cpp
+++ b/src/GksMeshAdapter/MeshCell.cpp
@@ -21,6 +21,8 @@ MeshCell::MeshCell(){
     isGhostCell = false;
 
     isWall = false;
+
+    isFluxBC = false;
 }
 
 bool MeshCell::isCoarseGhostCell()
diff --git a/src/GksMeshAdapter/MeshCell.h b/src/GksMeshAdapter/MeshCell.h
index 38d61836fe8d574dd59185648788b1752d9b41e1..ac0b9d40c53bfc365de9c0ea6738bb0551fbc9cc 100644
--- a/src/GksMeshAdapter/MeshCell.h
+++ b/src/GksMeshAdapter/MeshCell.h
@@ -60,6 +60,8 @@ struct VF_PUBLIC MeshCell{
 
     bool isWall;
 
+    bool isFluxBC;
+
     char type;
 
     MeshCell();
diff --git a/targets/apps/GKS/Flame7cm/Flame7cm.cpp b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
index d7cbe94015922651ce9b8f89479479eda1058966..598c1d6d134c8e11f1dc4488f6f389946a7d3e32 100644
--- a/targets/apps/GKS/Flame7cm/Flame7cm.cpp
+++ b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
@@ -4,6 +4,7 @@
 #include <math.h>
 #include <string>
 #include <iostream>
+#include <iomanip>
 #include <exception>
 #include <fstream>
 #include <memory>
@@ -42,6 +43,7 @@
 #include "GksGpu/BoundaryConditions/PassiveScalarDiriclet.h"
 #include "GksGpu/BoundaryConditions/InflowComplete.h"
 #include "GksGpu/BoundaryConditions/Open.h"
+#include "GksGpu/BoundaryConditions/Inflow.h"
 
 #include "GksGpu/Interface/Interface.h"
 #include "GksGpu/TimeStepping/NestedTimeStep.h"
@@ -187,14 +189,17 @@ 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<AdiabaticWall>( dataBase, Vec3(0, 0, 0), false );
+    //SPtr<BoundaryCondition> bcPX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), false );
+
+    SPtr<BoundaryCondition> bcMX_2 = std::make_shared<Inflow>( dataBase, Vec3( 100.0*U, 0.0, 0.0), prim.lambda, rho, 1.0, 0.0, 0.0 );
+    SPtr<BoundaryCondition> bcPX_2 = std::make_shared<Inflow>( dataBase, Vec3(-100.0*U, 0.0, 0.0), prim.lambda, rho, 1.0, 0.0, 0.0 );
 
     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; } );
+    bcMX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L && center.z < 0.25*H; } );
+    bcPX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x >  0.5*L && center.z < 0.25*H; } );
 
     //////////////////////////////////////////////////////////////////////////
     
@@ -227,19 +232,20 @@ void thermalCavity( std::string path, std::string simulationName )
 
     SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
     
-    bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0; } );
+    bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0 /*&& std::sqrt(center.x*center.x + center.y*center.y) >= 0.5*0.071*/; } );
     bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z > H  ; } );
 
     //////////////////////////////////////////////////////////////////////////
 
     //SPtr<BoundaryCondition> burner = std::make_shared<IsothermalWall>( dataBase, Vec3(0.0, 0.0, 0.0), 0.5*prim.lambda,  0.0, true );
 
-    SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, U, prim.lambda, 1.0, 0.0) );
+    SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, U, prim.lambda, 1.0, 1.0) );
     //SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, 10.0 * U, prim.lambda, 0.0, 0.0) );
 
     burner->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ 
 
-        return center.z < 0.0 && std::sqrt(center.x*center.x + center.y*center.y) < 0.5*0.071;
+        //return center.z < 0.0 && std::sqrt(center.x*center.x + center.y*center.y) < 0.5*0.071;
+        return center.z < 0.0 && std::sqrt(center.x*center.x) < 0.5*0.071;
     } );
 
     //////////////////////////////////////////////////////////////////////////
@@ -331,7 +337,7 @@ void thermalCavity( std::string path, std::string simulationName )
 
         if( 
             //( iter >= 2000 && iter % 100 == 0 ) || 
-            ( iter % 1000 == 0 )
+            ( iter % 10 == 0 )
           )
         {
             for( uint level = 0; level < dataBase->numberOfLevels; level++ )