diff --git a/CMakeLists.txt b/CMakeLists.txt
index d4a0c1a66c2ef816019426120c23df87a9175974..ca55cb0c5e3880200de5156cff317d087b19e539 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -160,7 +160,7 @@ IF (HULC.BUILD_VF_GKS)
     add_subdirectory(targets/apps/GKS/ConfinedCombustion)
     add_subdirectory(targets/apps/GKS/MethaneFlame)
     
-    #add_subdirectory(targets/apps/GKS/Room)
+    add_subdirectory(targets/apps/GKS/Room)
     add_subdirectory(targets/apps/GKS/RoomMultiGPU)
 ELSE()
   MESSAGE( STATUS "exclude Virtual Fluids GKS." )
diff --git a/src/GksGpu/CellUpdate/CellUpdate.cu b/src/GksGpu/CellUpdate/CellUpdate.cu
index f13f3abe14185eaa29cea3a42fda8ccb13227b25..9950f6dd5ed36a9ef10c670af0c1fbe79acbefc4 100644
--- a/src/GksGpu/CellUpdate/CellUpdate.cu
+++ b/src/GksGpu/CellUpdate/CellUpdate.cu
@@ -244,11 +244,11 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
         real Y_N2_ambient = 0.767;
         real Y_O2_ambient = 0.233;
 
-        real m_CH4 = 16.0;
-        real m_O2  = 32.0;
-        real m_N2  = 28.0;
-        real m_H2O = 34.0;
-        real m_CO2 = 44.0;
+        real m_CH4 = 16.0;  // g / mol
+        real m_O2  = 32.0;  // g / mol
+        real m_N2  = 28.0;  // g / mol
+        real m_H2O = 34.0;  // g / mol
+        real m_CO2 = 44.0;  // g / mol
 
         ///////////////////////////////////////////////////////////////////////////////
 
@@ -268,46 +268,24 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
 
         real Y_S;   // currently not modeled
 
-        //if( Z1 < 0.0 ) Z1 = 0.0;
-        //if( Z2 < 0.0 ) Z2 = 0.0;
-
         //////////////////////////////////////////////////////////////////////////
 
-        if( Y_CH4 > 1.0e-6 && Y_O2 > 1.0e-6 )
-
         {
-            const real heatOfReaction = 2.0e3;
+            const real heatOfReaction = 8.0e8;
 
             real s = m_CH4 / (2.0 * m_O2);      // refers to page 49 in FDS technical reference guide
 
-            real releasedHeat = rho * fminf(Y_CH4, s * Y_O2) * heatOfReaction;
-
-            if (Y_CH4 < s * Y_O2)
-            {
-                Y_CH4 = 0.0;
-                //Y_O2  = Y_O2  - Y_CH4 / s;
-            }
-            else
-            {
-                Y_CH4 = Y_CH4 - s * Y_O2;
-                //Y_O2  = 0.0;
-            }
+            real releasedHeat = rho * fminf(Y_CH4, s * Y_O2) / m_CH4 * heatOfReaction;
 
-            //if( Z2 < 0.0 ) printf( "%f, %f, %f \n", Z2, Y_O2_ambient, Y_O2 );
-
-            //Y_CH4 = Y_CH4 - fminf( Y_CH4    , s * Y_O2 );
-            //Y_O2  = Y_O2  - fminf( Y_CH4 / s,     Y_O2 );
+            if (Y_CH4 < s * Y_O2) Y_CH4 = 0.0;
+            else                  Y_CH4 = Y_CH4 - s * Y_O2;
 
             ///////////////////////////////////////////////////////////////////////////////
 
             real dZ1 = Z1 - Y_CH4 / Y_CH4_Inflow;
 
-            //if( dZ1 < 0.0 ) printf( "%f, %f, %f \n", Z1, Y_CH4, Y_CH4_Inflow );
-
             Z1 = Y_CH4 / Y_CH4_Inflow;
 
-            //Z2 = ( ( 1.0 - Z1 ) * Y_O2_ambient - Y_O2 ) / ( 1.0 + 2.0 * ( m_O2  / m_CH4 ) * Y_CH4_Inflow );
-
             Z2 = Z2 + dZ1;
 
             ///////////////////////////////////////////////////////////////////////////////
@@ -315,7 +293,7 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
             dataBase.data[RHO_S_1(cellIndex, dataBase.numberOfCells)] = Z1 * updatedConserved.rho;
             dataBase.data[RHO_S_2(cellIndex, dataBase.numberOfCells)] = Z2 * updatedConserved.rho;
 
-            //dataBase.data[ RHO_E(cellIndex, dataBase.numberOfCells) ] += releasedHeat;
+            dataBase.data[ RHO_E(cellIndex, dataBase.numberOfCells) ] = updatedConserved.rhoE + releasedHeat;
         }
     }
 
diff --git a/targets/apps/GKS/ConfinedCombustion/ConfinedCombustion.cpp b/targets/apps/GKS/ConfinedCombustion/ConfinedCombustion.cpp
index 2a16b420ee4ba0fe13b7ffe2125ce3cdc8014b0f..94a448913dece2679c708eeb206d0e9197c08fe8 100644
--- a/targets/apps/GKS/ConfinedCombustion/ConfinedCombustion.cpp
+++ b/targets/apps/GKS/ConfinedCombustion/ConfinedCombustion.cpp
@@ -60,7 +60,7 @@ void thermalCavity( std::string path, std::string simulationName )
     real Ba  = 0.1;
     real eps = 1.2;
     real Pr  = 0.71;
-    real K   = 8.0;
+    real K   = 3.0;
     
     real g   = 9.81;
     real rho = 1.2;
@@ -72,8 +72,8 @@ void thermalCavity( std::string path, std::string simulationName )
 				   + S_2               * 8.31445984848 / 32.00e-3      // CH4
 		           + (1.0 - S_1 - S_2) * 8.31445984848 / 28.00e-3;     // N2
 
-    real lambdaCold = 0.5 / ( R_Mixture *  300 );
-    real lambdaHot  = 0.5 / ( R_Mixture * 1200 );
+    real lambdaCold = 0.5 / ( 287 *  300 );
+    real lambdaHot  = 0.5 / ( 287 * 1200 );
     
     real mu = 1.0e-1;
 
@@ -82,7 +82,7 @@ void thermalCavity( std::string path, std::string simulationName )
 
     //real CFL = 0.8;
 
-    real dt  = 1.0e-4; //CFL * ( dx / ( ( U + cs ) * ( one + ( two * mu ) / ( U * dx * rho ) ) ) );
+    real dt  = 5.0e-6; //CFL * ( dx / ( ( U + cs ) * ( one + ( two * mu ) / ( U * dx * rho ) ) ) );
 
     *logging::out << logging::Logger::INFO_HIGH << "dt = " << dt << " s\n";
     //*logging::out << logging::Logger::INFO_HIGH << "U  = " << U  << " m/s\n";
@@ -159,6 +159,8 @@ void thermalCavity( std::string path, std::string simulationName )
     SPtr<BoundaryCondition> bcPX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0.0, 0.0, 0.0), false );
     //SPtr<BoundaryCondition> bcMX = std::make_shared<IsothermalWall>( dataBase, Vec3(0.0, 0.0, 0.0), lambdaCold, false );
     //SPtr<BoundaryCondition> bcPX = std::make_shared<IsothermalWall>( dataBase, Vec3(0.0, 0.0, 0.0), lambdaCold,  0.0, false );
+    //SPtr<BoundaryCondition> bcMX = std::make_shared<Pressure>( dataBase, 0.5 * rho / lambdaCold );
+    //SPtr<BoundaryCondition> bcPX = std::make_shared<Pressure>( dataBase, 0.5 * rho / lambdaCold );
 
     bcMX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L; } );
     bcPX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x >  0.5*L; } );
@@ -273,7 +275,7 @@ void thermalCavity( std::string path, std::string simulationName )
 
     cupsAnalyzer.start();
 
-    for( uint iter = 1; iter <= 100000000; iter++ )
+    for( uint iter = 1; iter <= 400000; iter++ )
     {
         //if( iter < 100000 )
         //{
@@ -294,7 +296,7 @@ void thermalCavity( std::string path, std::string simulationName )
             //( iter < 10000    && iter % 1000  == 0 ) ||
             //( iter < 1000000   && iter % 10000  == 0 ) ||
             //( iter < 10000000 && iter % 100000 == 0 )
-            ( iter < 20000 && iter % 100 == 0 )
+            ( iter <= 400000 && iter % 1000 == 0 )
           )
         {
             dataBase->copyDataDeviceToHost();