diff --git a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
index b86dc72e535c08860ef40b09de99b84c86da3a57..c605bb68e907ddaef16555bc87dfe6d2e68a014b 100644
--- a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
+++ b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
@@ -1,4 +1,4 @@
-#include "PointTimeseriesAnalyzer.h"
+#include "PointTimeSeriesAnalyzer.h"
 
 #include <cuda.h>
 #include <cuda_runtime.h>
diff --git a/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp b/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
index 0c2a8ece4a09489a9e1224b8d469841db20754ad..d1bc55c342db6f3b70f8a81fdff91022e59f6ec8 100644
--- a/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
+++ b/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
@@ -1,11 +1,11 @@
-#include "PointTimeseriesCollector.h"
+#include "PointTimeSeriesCollector.h"
 
 #include <iomanip>
 #include <fstream>
 
 #include "Core/Logger/Logger.h"
 
-#include "Analyzer/PointTimeseriesAnalyzer.h"
+#include "Analyzer/PointTimeSeriesAnalyzer.h"
 
 #include "Parameters/Parameters.h"
 
diff --git a/src/GksGpu/CellUpdate/Reaction.cuh b/src/GksGpu/CellUpdate/Reaction.cuh
index 87f548b7d0579e52409db70dae5fc56f56271582..5f7820d3137d7285e01eb5a981b0755bb2399159 100644
--- a/src/GksGpu/CellUpdate/Reaction.cuh
+++ b/src/GksGpu/CellUpdate/Reaction.cuh
@@ -19,6 +19,92 @@
 
 #include "CudaUtility/CudaRunKernel.hpp"
 
+inline __host__ __device__ real getTurbulentViscosityDeardorff(const DataBaseStruct& dataBase, const Parameters& parameters, const uint cellIndex, const ConservedVariables& cons )
+{
+    // See FDS 6 Technical Reference Guide, Section 4.2.3
+
+    PrimitiveVariables prim = toPrimitiveVariables(cons, parameters.K);
+
+    ConservedVariables neighborCons;
+    PrimitiveVariables neighborPrim;
+
+    real kSGS = c0o1;
+
+    {
+        real uHead = c1o2 * prim.U;
+
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 0, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            uHead += c1o4 * neighborPrim.U;
+        }
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 1, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            uHead += c1o4 * neighborPrim.U;
+        }
+
+        kSGS += c1o2 * ( prim.U - uHead ) * ( prim.U - uHead );
+    }
+
+    {
+        real vHead = c1o2 * prim.V;
+
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 2, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            vHead += c1o4 * neighborPrim.V;
+        }
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 3, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            vHead += c1o4 * neighborPrim.V;
+        }
+
+        kSGS += c1o2 * ( prim.V - vHead ) * ( prim.V - vHead );
+    }
+
+    {
+        real wHead = c1o2 * prim.W;
+
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 4, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            wHead += c1o4 * neighborPrim.W;
+        }
+        {
+            uint neighborCellIndex = dataBase.cellToCell[CELL_TO_CELL(cellIndex, 5, dataBase.numberOfCells)];
+            readCellData(cellIndex, dataBase, neighborCons);
+            neighborPrim = toPrimitiveVariables(neighborCons, parameters.K);
+
+            wHead += c1o4 * neighborPrim.W;
+        }
+
+        kSGS += c1o2 * ( prim.W - wHead ) * ( prim.W - wHead );
+    }
+
+    //real turbulentViscosity = prim.rho * parameters.dx * c1o10 * sqrt(kSGS) / 0.3;
+
+    dataBase.diffusivity[cellIndex] = (realAccumulator) kSGS;
+
+    //printf("%f", kSGS);
+
+    return kSGS;
+}
+
+
+
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 __host__ __device__ inline void chemicalReaction(DataBaseStruct dataBase, Parameters parameters, uint cellIndex, ConservedVariables& cons)
@@ -35,6 +121,8 @@ __host__ __device__ inline void chemicalReaction(DataBaseStruct dataBase, Parame
 
         //////////////////////////////////////////////////////////////////////////
 
+        //real diffusivity = getTurbulentViscosityDeardorff(dataBase, parameters, cellIndex, cons);
+        //real diffusivity = dataBase.diffusivity[ cellIndex ];
         real diffusivity = dataBase.diffusivity[ cellIndex ] / ( c6o1 * parameters.dx * parameters.dx * parameters.dt );
         dataBase.diffusivity[ cellIndex ] = c0o1;
 
@@ -42,10 +130,16 @@ __host__ __device__ inline void chemicalReaction(DataBaseStruct dataBase, Parame
 
         real mixingTimeScale = real(0.1) * parameters.dx * parameters.dx / diffusivity;
 
-        //real mixingTimeScale = parameters.dt;
+        //real kSGS = getTurbulentViscosityDeardorff(dataBase, parameters, cellIndex, cons);
 
-        //if( mixingTimeScale < one )
-        //    mixingTimeScale = one;
+        //real mixingTimeScale_d = parameters.dx * parameters.dx / parameters.D;
+
+        //real mixingTimeScale_u = real(0.4) * parameters.dx / sqrt( c2o3 * kSGS );
+
+        //real mixingTimeScale_g = sqrt( c2o1 * parameters.dx / fabs( parameters.force.z ) );
+
+        //real mixingTimeScale = fminf( mixingTimeScale_d, mixingTimeScale_u );
+        //mixingTimeScale      = fminf( mixingTimeScale_g, mixingTimeScale   );
 
         //////////////////////////////////////////////////////////////////////////
 
@@ -90,4 +184,4 @@ __host__ __device__ inline void chemicalReaction(DataBaseStruct dataBase, Parame
     }
 
 #endif // USE_PASSIVE_SCALAR
-}
+}
\ No newline at end of file
diff --git a/src/GksGpu/DataBase/DataBase.h b/src/GksGpu/DataBase/DataBase.h
index 847da7ea4cff773166674ec4368930a463539baa..e6b5dd6d33f503dbf32d561ed0cdf5d044e5778d 100644
--- a/src/GksGpu/DataBase/DataBase.h
+++ b/src/GksGpu/DataBase/DataBase.h
@@ -85,7 +85,7 @@ struct VF_PUBLIC DataBase : public std::enable_shared_from_this<DataBase>
     char* faceOrientation;
 
     uint* fineToCoarse;   // 9
-    uint* coarseToFine;   // 15
+    uint* coarseToFine;   // 9
 
     //////////////////////////////////////////////////////////////////////////
     // Host/Device data - READ MODIFY
diff --git a/src/GksGpu/Definitions/PassiveScalar.h b/src/GksGpu/Definitions/PassiveScalar.h
index b71ada48a0c05eaa726abae31bf76f7b56a2c0fd..97fca44e5285a9b9850ba5690eae0a420ce60449 100644
--- a/src/GksGpu/Definitions/PassiveScalar.h
+++ b/src/GksGpu/Definitions/PassiveScalar.h
@@ -1,6 +1,6 @@
 #ifndef PassiveScalar_H
 #define PassiveScalar_H
 
-//#define USE_PASSIVE_SCALAR
+#define USE_PASSIVE_SCALAR
 
 #endif
diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index d5dfcea381199f9aca5d4e08f43bc70af62b106e..087c6e9d2f4252402dd199afde275c139c52d1c5 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -284,6 +284,19 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
             real dTdx1 = dEdx1 - c2o1 * facePrim.U * dUdx1 - c2o1 * facePrim.V * dVdx1 - c2o1 * facePrim.W * dWdx1;
             real dTdx2 = dEdx2 - c2o1 * facePrim.U * dUdx2 - c2o1 * facePrim.V * dVdx2 - c2o1 * facePrim.W * dWdx2;
             real dTdx3 = dEdx3 - c2o1 * facePrim.U * dUdx3 - c2o1 * facePrim.V * dVdx3 - c2o1 * facePrim.W * dWdx3;
+    
+            //real E = c1o2 * ( facePrim.U * facePrim.U 
+            //                + facePrim.V * facePrim.V 
+            //                + facePrim.W * facePrim.W 
+            //                + ( parameters.K + c3o1 ) / ( c4o1 * facePrim.lambda ) );
+
+            //real dEdx1 = ( gradN.rhoE  - E * gradN.rho  );
+            //real dEdx2 = ( gradT1.rhoE - E * gradT1.rho );
+            //real dEdx3 = ( gradT2.rhoE - E * gradT2.rho );
+
+            //real dTdx1 = dEdx1 - facePrim.U * dUdx1 - facePrim.V * dVdx1 - facePrim.W * dWdx1;
+            //real dTdx2 = dEdx2 - facePrim.U * dUdx2 - facePrim.V * dVdx2 - facePrim.W * dWdx2;
+            //real dTdx3 = dEdx3 - facePrim.U * dUdx3 - facePrim.V * dVdx3 - facePrim.W * dWdx3;
 
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
             // this one works for some time
diff --git a/targets/apps/GKS/SalinasVazquez/SalinasVazquez.cpp b/targets/apps/GKS/SalinasVazquez/SalinasVazquez.cpp
index 80cd98ba5f26fe3688f4b7c9f5c9314b9bbcf64f..ecb96f7c90f080dace051e1c4d309073afec294d 100644
--- a/targets/apps/GKS/SalinasVazquez/SalinasVazquez.cpp
+++ b/targets/apps/GKS/SalinasVazquez/SalinasVazquez.cpp
@@ -100,7 +100,7 @@ void simulation( std::string path, std::string simulationName, uint restartIter
 
     real CFL = 0.25;
 
-    real dt  = CFL * ( dx / ( ( U + cs ) * ( one + ( two * mu ) / ( U * dx * rho ) ) ) );
+    real dt  = CFL * ( dx / ( ( U + cs ) * ( c1o1 + ( c2o1 * mu ) / ( U * dx * rho ) ) ) );
 
     *logging::out << logging::Logger::INFO_HIGH << "dt = " << dt << " s\n";
     *logging::out << logging::Logger::INFO_HIGH << "U  = " << U  << " s\n";
@@ -439,7 +439,7 @@ void simulation( std::string path, std::string simulationName, uint restartIter
             writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
         }
 
-        cupsAnalyzer.run( iter );
+        cupsAnalyzer.run( iter, dt );
 
         convergenceAnalyzer.run( iter );
 
@@ -545,15 +545,15 @@ int main( int argc, char* argv[])
     }
     catch (const std::exception& e)
     {     
-        *logging::out << logging::Logger::ERROR << e.what() << "\n";
+        *logging::out << logging::Logger::LOGGER_ERROR << e.what() << "\n";
     }
     catch (const std::bad_alloc& e)
     {  
-        *logging::out << logging::Logger::ERROR << "Bad Alloc:" << e.what() << "\n";
+        *logging::out << logging::Logger::LOGGER_ERROR << "Bad Alloc:" << e.what() << "\n";
     }
     catch (...)
     {
-        *logging::out << logging::Logger::ERROR << "Unknown exception!\n";
+        *logging::out << logging::Logger::LOGGER_ERROR << "Unknown exception!\n";
     }
 
     logFile.close();
diff --git a/targets/apps/GKS/SandiaFlame_1m/SandiaFlame_1m.cpp b/targets/apps/GKS/SandiaFlame_1m/SandiaFlame_1m.cpp
index 7a144be12d2d911cdd2a42133faf03de12221f0f..a57824fa89349b33ab2e166aa200957b4ca6110f 100644
--- a/targets/apps/GKS/SandiaFlame_1m/SandiaFlame_1m.cpp
+++ b/targets/apps/GKS/SandiaFlame_1m/SandiaFlame_1m.cpp
@@ -515,8 +515,8 @@ int main( int argc, char* argv[])
     //uint restartIter = 90000;
         
     uint gpuIndex = 1;
-    uint testIndex = 17;
-    uint nx = 128;
+    uint testIndex = 24;
+    uint nx = 64;
 
     if( argc > 1 ) gpuIndex    = atoi( argv[1] );
 
diff --git a/targets/apps/GKS/ThermalCavityMultiGPU/ThermalCavityMultiGPU.cpp b/targets/apps/GKS/ThermalCavityMultiGPU/ThermalCavityMultiGPU.cpp
index 23a4062d6b6a2cb568a9fa0e1cc4d08d1ca835aa..a72951916385d19736037d842e7df88c2c8fcb2e 100644
--- a/targets/apps/GKS/ThermalCavityMultiGPU/ThermalCavityMultiGPU.cpp
+++ b/targets/apps/GKS/ThermalCavityMultiGPU/ThermalCavityMultiGPU.cpp
@@ -51,7 +51,7 @@
 #include "GksGpu/Analyzer/CupsAnalyzer.h"
 #include "GksGpu/Analyzer/ConvergenceAnalyzer.h"
 #include "GksGpu/Analyzer/TurbulenceAnalyzer.h"
-#include "GksGpu/Analyzer/PointTimeseriesCollector.h"
+#include "GksGpu/Analyzer/PointTimeSeriesCollector.h"
 
 #include "GksGpu/Restart/Restart.h"
 
@@ -60,7 +60,7 @@
 //uint deviceMap [2] = {2,3};
 uint deviceMap [2] = {0,1};
 
-void simulation( std::string path, std::string simulationName, bool fine, uint restartIter )
+void simulation( std::string path, std::string simulationName, bool fine, bool highAspect, uint restartIter )
 {
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -77,7 +77,11 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     real L = 1.0;
-    real H = 0.25;
+    real W = 0.25;
+
+    real H = L;
+
+    if(highAspect) H = 2.0 * L;
 
     real dx = L / real(nx);
 
@@ -155,31 +159,30 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     if( mpiWorldSize == 2 )
     {
         startY = 0.0;
-        endY   = H;
+        endY   = W;
     }
     else
     {
-        startY =  rank / 2        * H - 3.0 * dx;
-        endY   = (rank / 2 + 1.0) * H + 3.0 * dx;
+        startY =  rank / 2        * W - 3.0 * dx;
+        endY   = (rank / 2 + 1.0) * W + 3.0 * dx;
     }
 
-    startZ = -0.5 * L;
-    endZ   =  0.5 * L;
+    startZ = - 0.5 * H;
+    endZ   =   0.5 * H;
 
     gridBuilder->addCoarseGrid(startX, startY, startZ,  
                                endX  , endY  , endZ  , dx);
 
+    std::cout << __LINE__ << std::endl;
+
     //////////////////////////////////////////////////////////////////////////
 
-    //real refL[4] = { 0.30, 0.45, 0.49, 0.4975  };
-    real refL[4] = { 0.30, 0.45, 0.475, 0.495  };
+    real refL[4] = { 0.2, 0.05, 0.025, 0.005 };
 
     if( fine )
     {
-        refL[1] = 0.4;
-        refL[2] = 0.45;
-        //refL[1] += 2;
-        //refL[2] += 2;
+        refL[1] = 0.1;
+        refL[2] = 0.05;
     }
 
     gridBuilder->setNumberOfLayers(6,6);
@@ -188,14 +191,14 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     Conglomerate coarseRefLevel;
 
-    if( rank % 2 == 0 ) coarseRefLevel.add( new Cuboid (-100.0,   -100.0, -100.0, 
-                                                        -refL[0],  100.0,  100.0 ) );
-    else                coarseRefLevel.add( new Cuboid ( refL[0], -100.0, -100.0, 
-                                                         100.0,    100.0,  100.0 ) );
+    if( rank % 2 == 0 ) coarseRefLevel.add( new Cuboid (-100.0,           -100.0, -100.0, 
+                                                        -0.5*L + refL[0],  100.0,  100.0 ) );
+    else                coarseRefLevel.add( new Cuboid ( 0.5*L - refL[0], -100.0, -100.0, 
+                                                         100.0,            100.0,  100.0 ) );
 
     coarseRefLevel.add( new Cuboid (-100.0, -100.0, -100.0,   
-                                     100.0,  100.0, -refL[0] ) );
-    coarseRefLevel.add( new Cuboid (-100.0, -100.0,  refL[0], 
+                                     100.0,  100.0, -0.5*H + refL[0] ) );
+    coarseRefLevel.add( new Cuboid (-100.0, -100.0,  0.5*H - refL[0], 
                                      100.0,  100.0,  100.0   ) );
 
     gridBuilder->addGrid( &coarseRefLevel, 1);
@@ -204,14 +207,14 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     Conglomerate firstRefLevel;
 
-    if( rank % 2 == 0 ) firstRefLevel.add( new Cuboid (-100.0,   -100.0, -100.0, 
-                                                       -refL[1],  100.0,  100.0 ) );
-    else                firstRefLevel.add( new Cuboid ( refL[1], -100.0, -100.0, 
-                                                        100.0,    100.0,  100.0 ) );
+    if( rank % 2 == 0 ) firstRefLevel.add( new Cuboid (-100.0,           -100.0, -100.0, 
+                                                       -0.5*L + refL[1],  100.0,  100.0 ) );
+    else                firstRefLevel.add( new Cuboid ( 0.5*L - refL[1], -100.0, -100.0, 
+                                                        100.0,            100.0,  100.0 ) );
 
     firstRefLevel.add( new Cuboid (-100.0, -100.0, -100.0,   
-                                    100.0,  100.0, -refL[1] ) );
-    firstRefLevel.add( new Cuboid (-100.0, -100.0,  refL[1], 
+                                    100.0,  100.0, -0.5*H + refL[1] ) );
+    firstRefLevel.add( new Cuboid (-100.0, -100.0,  0.5*H - refL[1], 
                                     100.0,  100.0,  100.0   ) );
 
     gridBuilder->addGrid( &firstRefLevel, 2);
@@ -220,20 +223,20 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     Conglomerate secondRefLevel;
 
-    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,   -100.0, -100.0, 
-                                                        -refL[2],  100.0,  100.0 ) );
-    else                secondRefLevel.add( new Cuboid ( refL[2], -100.0, -100.0, 
-                                                         100.0,    100.0,  100.0 ) );
+    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,           -100.0, -100.0, 
+                                                        -0.5*L + refL[2],  100.0,  100.0 ) );
+    else                secondRefLevel.add( new Cuboid ( 0.5*L - refL[2], -100.0, -100.0, 
+                                                         100.0,            100.0,  100.0 ) );
 
-    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,   -100.0, -100.0,   
-                                                        -refL[0],  100.0, -refL[2] ) );
-    else                secondRefLevel.add( new Cuboid ( refL[0], -100.0, -100.0,   
-                                                         100.0,    100.0, -refL[2] ) );
+    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,           -100.0, -100.0,   
+                                                        -0.5*L + refL[0],  100.0, -0.5*H + refL[2] ) );
+    else                secondRefLevel.add( new Cuboid ( 0.5*L - refL[0], -100.0, -100.0,   
+                                                         100.0,            100.0, -0.5*H + refL[2] ) );
 
-    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,   -100.0,  refL[2], 
-                                                        -refL[0],  100.0,  100.0   ) );
-    else                secondRefLevel.add( new Cuboid ( refL[0], -100.0,  refL[2], 
-                                                         100.0,    100.0,  100.0   ) );
+    if( rank % 2 == 0 ) secondRefLevel.add( new Cuboid (-100.0,           -100.0,  0.5*H - refL[2], 
+                                                        -0.5*L + refL[0],  100.0,  100.0   ) );
+    else                secondRefLevel.add( new Cuboid ( 0.5*L - refL[0], -100.0,  0.5*H - refL[2], 
+                                                         100.0,            100.0,  100.0   ) );
 
     gridBuilder->addGrid( &secondRefLevel, 3);
 
@@ -241,10 +244,10 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     Conglomerate thirdRefLevel;
 
-    if( rank % 2 == 0 ) thirdRefLevel.add( new Cuboid (-100.0,   -100.0, -100.0, 
-                                                        -refL[3],  100.0,  100.0 ) );
-    else                thirdRefLevel.add( new Cuboid ( refL[3], -100.0, -100.0, 
-                                                        100.0,    100.0,  100.0 ) );
+    if( rank % 2 == 0 ) thirdRefLevel.add( new Cuboid (-100.0,           -100.0, -100.0, 
+                                                       -0.5*L + refL[3],  100.0,  100.0 ) );
+    else                thirdRefLevel.add( new Cuboid ( 0.5*L - refL[3], -100.0, -100.0, 
+                                                        100.0,            100.0,  100.0 ) );
 
     if( fine ) gridBuilder->addGrid( &thirdRefLevel, 4);
 
@@ -262,8 +265,8 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     }
     else
     {
-        startY =   real(rank/2)         * H;
-        endY   = ( real(rank/2) + 1.0 ) * H;
+        startY =   real(rank/2)         * W;
+        endY   = ( real(rank/2) + 1.0 ) * W;
     }
 
     startZ = -100.0;
@@ -334,8 +337,8 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     SPtr<BoundaryCondition> bcMZ = std::make_shared<AdiabaticWall>( dataBase, Vec3(0,0,0), true );
     SPtr<BoundaryCondition> bcPZ = std::make_shared<AdiabaticWall>( dataBase, Vec3(0,0,0), true );
 
-    bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < -0.5*L; } );
-    bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z >  0.5*L; } );
+    bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < -0.5*H; } );
+    bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z >  0.5*H; } );
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -345,7 +348,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
         SPtr<BoundaryCondition> bcPY = std::make_shared<Periodic>(dataBase);
 
         bcMY->findBoundaryCells(meshAdapter, false, [&](Vec3 center) { return center.y < 0; });
-        bcPY->findBoundaryCells(meshAdapter, false, [&](Vec3 center) { return center.y > H; });
+        bcPY->findBoundaryCells(meshAdapter, false, [&](Vec3 center) { return center.y > W; });
 
         dataBase->boundaryConditions.push_back(bcMY);
         dataBase->boundaryConditions.push_back(bcPY);
@@ -377,10 +380,10 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     {
         Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> ConservedVariables {
 
-            real Th = 1.0 / lambdaHot;
-            real Tc = 1.0 / lambdaCold;
-            real T = Th - (Th - Tc)*((cellCenter.x + 0.5 * L) / L);
-            real lambdaLocal = 1.0 / T;
+            //real Th = 1.0 / lambdaHot;
+            //real Tc = 1.0 / lambdaCold;
+            //real T = Th - (Th - Tc)*((cellCenter.x + 0.5 * L) / L);
+            //real lambdaLocal = 1.0 / T;
 
             return toConservedVariables(PrimitiveVariables(rho, 0.0, 0.0, 0.0, lambda), parameters.K);
         });
@@ -431,17 +434,17 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     auto pointTimeSeriesCollector = std::make_shared<PointTimeSeriesCollector>();
 
-    for( real y = 0.5 * H; y < real( mpiWorldSize / 2 ) * H; y += H )
+    for( real y = 0.5 * W; y < real( mpiWorldSize / 2 ) * W; y += W )
     {
-        if( subDomainBox->isInside( -0.485, y, -0.3*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y, -0.3*L ), 'W', 10000 );
-        if( subDomainBox->isInside( -0.485, y, -0.1*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y, -0.1*L ), 'W', 10000 );
-        if( subDomainBox->isInside( -0.485, y,  0.1*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y,  0.1*L ), 'W', 10000 );
-        if( subDomainBox->isInside( -0.485, y,  0.3*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y,  0.3*L ), 'W', 10000 );
+        if( subDomainBox->isInside( -0.485, y, -0.3*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y, -0.3*H ), 'W', 10000 );
+        if( subDomainBox->isInside( -0.485, y, -0.1*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y, -0.1*H ), 'W', 10000 );
+        if( subDomainBox->isInside( -0.485, y,  0.1*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y,  0.1*H ), 'W', 10000 );
+        if( subDomainBox->isInside( -0.485, y,  0.3*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3( -0.485, y,  0.3*H ), 'W', 10000 );
         
-        if( subDomainBox->isInside(  0.485, y, -0.3*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y, -0.3*L ), 'W', 10000 );
-        if( subDomainBox->isInside(  0.485, y, -0.1*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y, -0.1*L ), 'W', 10000 );
-        if( subDomainBox->isInside(  0.485, y,  0.1*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y,  0.1*L ), 'W', 10000 );
-        if( subDomainBox->isInside(  0.485, y,  0.3*L ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y,  0.3*L ), 'W', 10000 );
+        if( subDomainBox->isInside(  0.485, y, -0.3*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y, -0.3*H ), 'W', 10000 );
+        if( subDomainBox->isInside(  0.485, y, -0.1*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y, -0.1*H ), 'W', 10000 );
+        if( subDomainBox->isInside(  0.485, y,  0.1*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y,  0.1*H ), 'W', 10000 );
+        if( subDomainBox->isInside(  0.485, y,  0.3*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y,  0.3*H ), 'W', 10000 );
     }
 
     //////////////////////////////////////////////////////////////////////////
@@ -503,6 +506,12 @@ int main( int argc, char* argv[])
 {
     //////////////////////////////////////////////////////////////////////////
 
+    bool fine = false;
+
+    bool highAspect = true;
+
+    //////////////////////////////////////////////////////////////////////////
+
 #ifdef _WIN32
     MPI_Init(&argc, &argv);
     int rank = 0;
@@ -528,7 +537,10 @@ int main( int argc, char* argv[])
     std::string path( "out/" );
 #endif
 
-    std::string simulationName ( "ThermalCavity3D_coarse" );
+    std::string simulationName ( "ThermalCavity3D" );
+
+    if(fine) simulationName += "_fine";
+    else     simulationName += "_coarse";
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -573,7 +585,7 @@ int main( int argc, char* argv[])
 
         if( argc > 1 ) restartIter = atoi( argv[1] );
 
-        simulation(path, simulationName, false, restartIter);
+        simulation(path, simulationName, fine, highAspect, restartIter);
     }
     catch (const std::exception& e)
     {