diff --git a/src/GksGpu/CellUpdate/CellUpdate.cu b/src/GksGpu/CellUpdate/CellUpdate.cu
index 97664ac1a3b35f4650c945da35eccb9a9824dbae..4e1cb8112607fc79493b27eb1e4085c94444796b 100644
--- a/src/GksGpu/CellUpdate/CellUpdate.cu
+++ b/src/GksGpu/CellUpdate/CellUpdate.cu
@@ -94,7 +94,7 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    if(false)
+    if(parameters.forcingSchemeIdx == 0)
     {
         // consistent source term treatment of Tian et al. (2007)
         cons.rhoU += parameters.force.x * parameters.dt * cons.rho;
@@ -109,7 +109,7 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
         dataBase.massFlux[VEC_Z(cellIndex, dataBase.numberOfCells)] = zero;
     }
 
-    if(false)
+    if(parameters.forcingSchemeIdx == 1)
     {
         // forcing only on density variation
         cons.rhoU += parameters.force.x * parameters.dt * ( cons.rho - parameters.rhoRef );
@@ -124,7 +124,7 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
         dataBase.massFlux[VEC_Z(cellIndex, dataBase.numberOfCells)] = zero;
     }
 
-    if(true)
+    if(parameters.forcingSchemeIdx == 2)
     {
         PrimitiveVariables prim = toPrimitiveVariables(cons, parameters.K);
         real lambda = prim.lambda;
diff --git a/src/GksGpu/CellUpdate/Reaction.cuh b/src/GksGpu/CellUpdate/Reaction.cuh
index 30bedd929b5289f70318b604596666852a4838fe..df39a4b19762156b83cbb930094a934fdb90eb8a 100644
--- a/src/GksGpu/CellUpdate/Reaction.cuh
+++ b/src/GksGpu/CellUpdate/Reaction.cuh
@@ -29,7 +29,7 @@ __host__ __device__ inline void chemicalReaction(DataBaseStruct dataBase, Parame
     {
         CellProperties cellProperties = dataBase.cellProperties[ cellIndex ];
 
-        if( !isCellProperties( cellProperties, CELL_PROPERTIES_FINE_GHOST ) )
+        //if( !isCellProperties( cellProperties, CELL_PROPERTIES_FINE_GHOST ) )
         //if( !isCellProperties( cellProperties, CELL_PROPERTIES_GHOST ) )
         {
             PrimitiveVariables prim = toPrimitiveVariables(cons, parameters.K);
diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index 1f20f06e03226fdca5f574f9757cfde944df83bd..a397b5af5fb47a79b111f04eb262483de5fe5619 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -116,25 +116,30 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
 
     //////////////////////////////////////////////////////////////////////////
 
-    //{   // SpongeLayer
-    //    real x = dataBase.faceCenter[ VEC_X(faceIndex, dataBase.numberOfFaces) ];
-    //    real z = dataBase.faceCenter[ VEC_Z(faceIndex, dataBase.numberOfFaces) ];
+    if( parameters.useSpongeLayer )
+    {
+        if( parameters.spongeLayerIdx == 0 )
+        {
+            real x = dataBase.faceCenter[VEC_X(faceIndex, dataBase.numberOfFaces)];
+            real z = dataBase.faceCenter[VEC_Z(faceIndex, dataBase.numberOfFaces)];
 
-    //    real muNew = parameters.mu;
+            real muNew = parameters.mu;
 
-    //    if( fabsf(x) > three )
-    //    {
-    //        muNew += ( fabsf(x) - three ) * ten * parameters.mu;
-    //    }
-    //    if( fabsf(z) > seven )
-    //    {
-    //        muNew += ( fabs(z) - seven ) * ten * parameters.mu;
-    //    }
+            //if (fabsf(x) > three)
+            //{
+            //    muNew += (fabsf(x) - three) * ten * parameters.mu;
+            //}
 
-    //    //parameters.Pr = muNew / parameters.mu;
+            real zStart = real(0.35);
 
-    //    parameters.mu = muNew;
-    //}
+            if (fabsf(z) > zStart)
+            {
+                muNew += (fabs(z) - zStart) * ten * ten * ten * parameters.mu;
+            }
+
+            parameters.mu = muNew;
+        }
+    }
 
     //////////////////////////////////////////////////////////////////////////
 
diff --git a/src/GksGpu/Parameters/Parameters.h b/src/GksGpu/Parameters/Parameters.h
index 4545ea3f423b2c6208cd35236d317eb8a592c40e..db640b877d2bd605392a218852d01689146c3118 100644
--- a/src/GksGpu/Parameters/Parameters.h
+++ b/src/GksGpu/Parameters/Parameters.h
@@ -40,6 +40,15 @@ struct  VF_PUBLIC Parameters
 
     //////////////////////////////////////////////////////////////////////////
 
+    bool useSpongeLayer = false;
+    uint spongeLayerIdx = 0;
+
+    //////////////////////////////////////////////////////////////////////////
+
+    uint forcingSchemeIdx = 0;
+
+    //////////////////////////////////////////////////////////////////////////
+
     bool enableReaction = false;
 
     real heatOfReaction = real(8000.0); // kJ / kmol  
diff --git a/src/GksMeshAdapter/MeshCell.h b/src/GksMeshAdapter/MeshCell.h
index 48f59cfc814a1a18de413378ae36abe7d81d0c5a..969f348c06a94ff23d52dfb59d8118b06e232418 100644
--- a/src/GksMeshAdapter/MeshCell.h
+++ b/src/GksMeshAdapter/MeshCell.h
@@ -55,8 +55,8 @@ struct VF_PUBLIC MeshCell{
 
     // order: +x, -x, +y, -y, +z, -z (see cellToCell)
     bool_6 faceExists;
-
-    bool isGhostCell;
+    
+    bool isGhostCell;   // this denotes cells that do not have all neighbors, excluding coarse ghost cells
 
     bool isWall;
 
diff --git a/targets/apps/GKS/Flame7cm/Flame7cm.cpp b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
index c4f4958b7eb3982d9701467c6205061ec450c58b..aa6cdd113c3e48d507a4e8d61e9933a858449383 100644
--- a/targets/apps/GKS/Flame7cm/Flame7cm.cpp
+++ b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
@@ -94,7 +94,8 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     real dt  = CFL * ( dx / ( ( U + cs ) * ( one + ( two * mu ) / ( U * dx * rho ) ) ) );
 
-    real dh = 4192.0; // kJ / kmol  / T_FAKTOR
+    //real dh = 4192.0; // kJ / kmol  / T_FAKTOR
+    real dh = 8000.0; // kJ / kmol  / T_FAKTOR
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -104,7 +105,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
     *logging::out << logging::Logger::INFO_HIGH << "mu = " << mu << " kg/sm\n";
     *logging::out << logging::Logger::INFO_HIGH << "Pr = " << Pr << "\n";
 
-    *logging::out << logging::Logger::INFO_HIGH << "HRR = " << U * rho * M_PI * R * R * ( dh * 100 ) / 0.016 / 1000.0 << " kW\n";
+    *logging::out << logging::Logger::INFO_HIGH << "HRR = " << U * /*rho*/ 0.68 * M_PI * R * R * ( dh * 100 ) / 0.016 / 1000.0 << " kW\n";
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -139,6 +140,13 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
     parameters.usePassiveScalarLimiter = true;
     parameters.useSmagorinsky          = false;
 
+    parameters.reactionLimiter = 1.005;
+
+    parameters.useSpongeLayer = true;
+    parameters.spongeLayerIdx = 0;
+
+    parameters.forcingSchemeIdx = 2;
+
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     auto gridFactory = GridFactory::make();
@@ -267,7 +275,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     //////////////////////////////////////////////////////////////////////////
 
-    SPtr<BoundaryCondition> burner = std::make_shared<CreepingMassFlux>( dataBase, rho, U, prim.lambda );
+    SPtr<BoundaryCondition> burner = std::make_shared<CreepingMassFlux>( dataBase, /*rho*/0.68, U, prim.lambda );
 
     burner->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ 
         
@@ -342,7 +350,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    CupsAnalyzer cupsAnalyzer( dataBase, true, 30.0 );
+    CupsAnalyzer cupsAnalyzer( dataBase, true, 30.0, true, 10000 );
 
     ConvergenceAnalyzer convergenceAnalyzer( dataBase, 10000 );
 
@@ -361,7 +369,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
         TimeStepping::nestedTimeStep(dataBase, parameters, 0);
 
         if( 
-            //( iter >= 100 && iter % 10 == 0 ) || 
+            //( iter >= 144820 && iter % 1 == 0 ) || 
             ( iter % 10000 == 0 )
           )
         {
@@ -423,7 +431,7 @@ int main( int argc, char* argv[])
     try
     {
         uint restartIter = INVALID_INDEX;
-        //uint restartIter = 100000;
+        //uint restartIter = 140000;
 
         if( argc > 1 ) restartIter = atoi( argv[1] );