From c952bf6fa2171d5afd35198b62efe37d7f601f54 Mon Sep 17 00:00:00 2001
From: "LEGOLAS\\lenz" <lenz@irmb.tu-bs.de>
Date: Fri, 27 Mar 2020 09:21:29 +0100
Subject: [PATCH] implements alternative SutherlandsLaw

---
 CMakeLists.txt                                |  27 +--
 src/GksGpu/Analyzer/HeatFluxAnalyzer.cu       |   8 +-
 src/GksGpu/Analyzer/KineticEnergyAnalyzer.cu  |   2 +
 ...Collector.h => PointTimeSeriesCollector.h} |   0
 .../BoundaryConditions/InflowComplete.cu      | 184 +++++++++---------
 src/GksGpu/FluxComputation/AssembleFlux.cuh   |  20 +-
 src/GksGpu/FluxComputation/SutherlandsLaw.cuh |  25 ++-
 src/GksGpu/Parameters/Parameters.h            |   3 +-
 .../RayleighBenardMultiGPU.cpp                |  80 ++++----
 9 files changed, 181 insertions(+), 168 deletions(-)
 rename src/GksGpu/Analyzer/{PointTimeseriesCollector.h => PointTimeSeriesCollector.h} (100%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6842086ba..d03261516 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -77,8 +77,8 @@ ENDIF(MSVC)
 ###                         OPTIONS                       ###
 #############################################################
 option(BUILD_SHARED_LIBS        "Build shared libraries"      ON )
-option(VF.BUILD_VF_GPU          "Build VirtualFluids GPU"     ON )
-option(VF.BUILD_VF_GKS          "Build VirtualFluids GKS"     OFF)
+option(VF.BUILD_VF_GPU          "Build VirtualFluids GPU"     OFF)
+option(VF.BUILD_VF_GKS          "Build VirtualFluids GKS"     ON )
 option(VF.BUILD_VF_TRAFFIC      "Build VirtualFluids Traffic" OFF)
 option(VF.BUILD_JSONCPP         "Builds json cpp "            OFF)
 option(VF.BUILD_NUMERIC_TESTS   "Build numeric tests"         OFF)
@@ -166,33 +166,34 @@ IF (VF.BUILD_VF_GKS)
     #add_subdirectory(targets/apps/GKS/LiFuXu)
 
 	#add_subdirectory(targets/apps/GKS/TaylorGreen3D)
-    add_subdirectory(targets/apps/GKS/DrivenCavity3D)
+    #add_subdirectory(targets/apps/GKS/DrivenCavity3D)
     #add_subdirectory(targets/apps/GKS/ThermalCavity)
 
-    add_subdirectory(targets/apps/GKS/ThermalCavityMultiGPU)
+    #add_subdirectory(targets/apps/GKS/ThermalCavityMultiGPU)
     #add_subdirectory(targets/apps/GKS/DrivenCavityMultiGPU)
+    add_subdirectory(targets/apps/GKS/RayleighBenardMultiGPU)
 
-    add_subdirectory(targets/apps/GKS/SalinasVazquez)
+    #add_subdirectory(targets/apps/GKS/SalinasVazquez)
     #add_subdirectory(targets/apps/GKS/BoundaryJet)
 
     #add_subdirectory(targets/apps/GKS/PropaneFlame)
-    add_subdirectory(targets/apps/GKS/ConfinedCombustion)
+    #add_subdirectory(targets/apps/GKS/ConfinedCombustion)
     #add_subdirectory(targets/apps/GKS/MethaneFlame)
     
     #add_subdirectory(targets/apps/GKS/Room)
     #add_subdirectory(targets/apps/GKS/RoomMultiGPU)
     #add_subdirectory(targets/apps/GKS/RoomFire)
-    add_subdirectory(targets/apps/GKS/RoomFireExtended)
-    add_subdirectory(targets/apps/GKS/ConcreteHeatFluxBCTest)
+    #add_subdirectory(targets/apps/GKS/RoomFireExtended)
+    #add_subdirectory(targets/apps/GKS/ConcreteHeatFluxBCTest)
     
     #add_subdirectory(targets/apps/GKS/PoolFire)
-    add_subdirectory(targets/apps/GKS/Flame7cm)
-    add_subdirectory(targets/apps/GKS/SandiaFlame_1m)
+    #add_subdirectory(targets/apps/GKS/Flame7cm)
+    #add_subdirectory(targets/apps/GKS/SandiaFlame_1m)
     #add_subdirectory(targets/apps/GKS/Candle)
     
-    add_subdirectory(targets/apps/GKS/MultiGPU)
-    add_subdirectory(targets/apps/GKS/MultiGPU_nD)
-    add_subdirectory(targets/apps/GKS/SingleGPU)
+    #add_subdirectory(targets/apps/GKS/MultiGPU)
+    #add_subdirectory(targets/apps/GKS/MultiGPU_nD)
+    #add_subdirectory(targets/apps/GKS/SingleGPU)
 ELSE()
   MESSAGE( STATUS "exclude Virtual Fluids GKS." )
 ENDIF()
diff --git a/src/GksGpu/Analyzer/HeatFluxAnalyzer.cu b/src/GksGpu/Analyzer/HeatFluxAnalyzer.cu
index e03c8fe52..266d3b564 100644
--- a/src/GksGpu/Analyzer/HeatFluxAnalyzer.cu
+++ b/src/GksGpu/Analyzer/HeatFluxAnalyzer.cu
@@ -106,13 +106,7 @@ __host__ __device__ void heatFluxFunction(DataBaseStruct& dataBase, GksGpu::Boun
 
     real r   = parameters.lambdaRef / lambda;
 
-    real mu;
-    if ( parameters.viscosityModel == ViscosityModel::constant ){
-        mu = parameters.mu;
-    }
-    else if ( parameters.viscosityModel == ViscosityModel::sutherlandsLaw ){
-        mu = sutherlandsLaw( parameters, r );
-    }
+    real mu = getViscosity(parameters, r);
 
     heatFlux[ startIndex + index ] = c1o4 * (parameters.K + c5o1) * ( mu / parameters.Pr ) / parameters.dx * ( c1o1 / domainPrim.lambda - c1o1 / ghostPrim.lambda );
 }
diff --git a/src/GksGpu/Analyzer/KineticEnergyAnalyzer.cu b/src/GksGpu/Analyzer/KineticEnergyAnalyzer.cu
index e2df78f78..70b130e6a 100644
--- a/src/GksGpu/Analyzer/KineticEnergyAnalyzer.cu
+++ b/src/GksGpu/Analyzer/KineticEnergyAnalyzer.cu
@@ -52,6 +52,8 @@ bool KineticEnergyAnalyzer::run(uint iter)
     this->kineticEnergyTimeSeries.push_back( EKin );
 
     //*logging::out << logging::Logger::INFO_HIGH << "EKin = " << EKin << "\n";
+
+    return true;
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/GksGpu/Analyzer/PointTimeseriesCollector.h b/src/GksGpu/Analyzer/PointTimeSeriesCollector.h
similarity index 100%
rename from src/GksGpu/Analyzer/PointTimeseriesCollector.h
rename to src/GksGpu/Analyzer/PointTimeSeriesCollector.h
diff --git a/src/GksGpu/BoundaryConditions/InflowComplete.cu b/src/GksGpu/BoundaryConditions/InflowComplete.cu
index b8e62893c..77ef127f3 100644
--- a/src/GksGpu/BoundaryConditions/InflowComplete.cu
+++ b/src/GksGpu/BoundaryConditions/InflowComplete.cu
@@ -202,98 +202,98 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    if( false )
-    {
-        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] = c0o1; 
-            ay[i] = c0o1; 
-            az[i] = c0o1; 
-            at[i] = c0o1;
-        }
-        
-        {
-            ConservedVariables gradN, gradT1, gradT2;
-
-            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);
-            applyFluxToNegCell(dataBase, ghostCellIdx , flux, 'z', parameters);
-        }
-    }
+    //if( false )
+    //{
+    //    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] = c0o1; 
+    //        ay[i] = c0o1; 
+    //        az[i] = c0o1; 
+    //        at[i] = c0o1;
+    //    }
+    //    
+    //    {
+    //        ConservedVariables gradN, gradT1, gradT2;
+
+    //        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);
+    //        applyFluxToNegCell(dataBase, ghostCellIdx , flux, 'z', parameters);
+    //    }
+    //}
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 }
 
diff --git a/src/GksGpu/FluxComputation/AssembleFlux.cuh b/src/GksGpu/FluxComputation/AssembleFlux.cuh
index 4b2773275..bda2ec6f0 100644
--- a/src/GksGpu/FluxComputation/AssembleFlux.cuh
+++ b/src/GksGpu/FluxComputation/AssembleFlux.cuh
@@ -205,19 +205,13 @@ __host__ __device__ inline void computeTimeDerivative( const PrimitiveVariables&
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-__host__ __device__ inline void computeTimeCoefficients(const PrimitiveVariables & facePrim, const Parameters& parameters, real timeCoefficients[4])
+__host__ __device__ inline void computeTimeCoefficients(const PrimitiveVariables & facePrim, Parameters& parameters, real timeCoefficients[4])
 {
     real r   = parameters.lambdaRef / facePrim.lambda;
 
     //if( r < zero ) printf( "ERROR: %f/%f\n", parameters.lambdaRef, facePrim.lambda );
 
-    real mu;
-    if ( parameters.viscosityModel == ViscosityModel::constant ){
-        mu = parameters.mu;
-    }
-    else if ( parameters.viscosityModel == ViscosityModel::sutherlandsLaw ){
-        mu = sutherlandsLaw( parameters, r );
-    }
+    real mu = getViscosity(parameters, r);
 
     real tau = c2o1 * facePrim.lambda * mu / facePrim.rho;
 
@@ -230,17 +224,11 @@ __host__ __device__ inline void computeTimeCoefficients(const PrimitiveVariables
     timeCoefficients[3] =                   tau     ;
 }
 
-__host__ __device__ inline void getTau(const PrimitiveVariables & facePrim, const Parameters& parameters, real& tau)
+__host__ __device__ inline void getTau(const PrimitiveVariables & facePrim, Parameters& parameters, real& tau)
 {
     real r   = parameters.lambdaRef / facePrim.lambda;
 
-    real mu;
-    if ( parameters.viscosityModel == ViscosityModel::constant ){
-        mu = parameters.mu;
-    }
-    else if ( parameters.viscosityModel == ViscosityModel::sutherlandsLaw ){
-        mu = sutherlandsLaw( parameters, r );
-    }
+    real mu = getViscosity(parameters, r);  mu = sutherlandsLaw2( parameters, r );
 
     tau = c2o1 * facePrim.lambda * mu / facePrim.rho;
 }
diff --git a/src/GksGpu/FluxComputation/SutherlandsLaw.cuh b/src/GksGpu/FluxComputation/SutherlandsLaw.cuh
index 4f1f94d2a..6373fee2d 100644
--- a/src/GksGpu/FluxComputation/SutherlandsLaw.cuh
+++ b/src/GksGpu/FluxComputation/SutherlandsLaw.cuh
@@ -15,7 +15,7 @@
 
 namespace GksGpu {
 
-inline __host__ __device__ real sutherlandsLaw(const Parameters & parameters, const real r)
+inline __host__ __device__ real sutherlandsLaw(Parameters & parameters, const real r)
 {
     real S  = real( 110.5 );
 
@@ -26,6 +26,29 @@ inline __host__ __device__ real sutherlandsLaw(const Parameters & parameters, co
     return parameters.mu * sqrt( r * r * r ) * ( C  + c1o1 ) / ( r  + C );
 }
 
+inline __host__ __device__ real sutherlandsLaw2(Parameters & parameters, const real r)
+{
+    real Smu = real( 0.648 );
+
+    real Sk  = real( 0.368 );
+
+    parameters.Pr *= ( ( Smu  + c1o1 ) / ( Sk  + c1o1 ) ) * ( ( r  + Sk ) / ( r  + Smu ) );
+
+    return parameters.mu * sqrt( r * r * r ) * ( Smu  + c1o1 ) / ( r  + Smu );
+}
+
+inline __host__ __device__ real getViscosity(Parameters & parameters, const real r)
+{
+    if ( parameters.viscosityModel == ViscosityModel::sutherlandsLaw ){
+        return sutherlandsLaw( parameters, r );
+    }
+    else if ( parameters.viscosityModel == ViscosityModel::sutherlandsLaw2 ){
+        return sutherlandsLaw2( parameters, r );
+    }
+
+    return parameters.mu;
+}
+
 } // namespace GksGpu
 
 #endif
diff --git a/src/GksGpu/Parameters/Parameters.h b/src/GksGpu/Parameters/Parameters.h
index c65023c20..ec1dae008 100644
--- a/src/GksGpu/Parameters/Parameters.h
+++ b/src/GksGpu/Parameters/Parameters.h
@@ -10,7 +10,8 @@ namespace GksGpu {
 
 enum class VF_PUBLIC ViscosityModel{
     constant,
-    sutherlandsLaw
+    sutherlandsLaw,
+    sutherlandsLaw2
 };
 
 struct  VF_PUBLIC Parameters
diff --git a/targets/apps/GKS/RayleighBenardMultiGPU/RayleighBenardMultiGPU.cpp b/targets/apps/GKS/RayleighBenardMultiGPU/RayleighBenardMultiGPU.cpp
index 333917e1d..03c14edda 100644
--- a/targets/apps/GKS/RayleighBenardMultiGPU/RayleighBenardMultiGPU.cpp
+++ b/targets/apps/GKS/RayleighBenardMultiGPU/RayleighBenardMultiGPU.cpp
@@ -99,8 +99,8 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     real dx = L / real(nx);
 
-    //real Ra = 1.0e7;
-    real Ra = 1.0e2;
+    real Ra = 3.0e6;
+    //real Ra = 1.0e2;
 
     real Ba  = 0.1;
     real eps = 0.8;
@@ -119,7 +119,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     real cs  = sqrt( ( ( K + 4.0 ) / ( K + 2.0 ) ) / ( 2.0 * lambda ) );
     real U   = sqrt( Ra ) * mu / ( rho * L );
 
-    real CFL = 0.05;
+    real CFL = 0.5;
 
     real dt  = CFL * ( dx / ( ( U + cs ) * ( c1o1 + ( c2o1 * mu ) / ( U * dx * rho ) ) ) );
 
@@ -129,7 +129,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     //////////////////////////////////////////////////////////////////////////
 
-    Parameters parameters;
+    GksGpu::Parameters parameters;
 
     parameters.K  = K;
     parameters.Pr = Pr;
@@ -144,8 +144,8 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     parameters.lambdaRef = lambda;
 
-    //parameters.viscosityModel = ViscosityModel::sutherlandsLaw;
-    parameters.viscosityModel = ViscosityModel::constant;
+    parameters.viscosityModel = GksGpu::ViscosityModel::sutherlandsLaw2;
+    //parameters.viscosityModel = ViscosityModel::constant;
 
     parameters.forcingSchemeIdx = 0;
 
@@ -237,7 +237,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     if( sideLengthZ == 1 || rankZ == 1 ) firstRefLevel.add( new Cuboid (-100.0, -100.0,  0.5*L - refL[1], 
                                                                          100.0,  100.0,  100.0           ) );
 
-    //gridBuilder->addGrid( &firstRefLevel, 2);
+    gridBuilder->addGrid( &firstRefLevel, 2);
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -346,7 +346,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    auto dataBase = std::make_shared<DataBase>( "GPU" );
+    auto dataBase = std::make_shared<GksGpu::DataBase>( "GPU" );
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -397,18 +397,18 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     dataBase->setCommunicators( meshAdapter );
 
-    CudaUtility::printCudaMemoryUsage();
+    GksGpu::CudaUtility::printCudaMemoryUsage();
 
     if( restartIter == INVALID_INDEX )
     {
-        Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> ConservedVariables {
+        GksGpu::Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> GksGpu::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;
 
-            return toConservedVariables(PrimitiveVariables(rho, 0.0, 0.0, 0.0, lambda), parameters.K);
+            return GksGpu::toConservedVariables(GksGpu::PrimitiveVariables(rho, 0.0, 0.0, 0.0, lambda), parameters.K);
         });
 
         if (rank == 0) writeVtkXMLParallelSummaryFile(dataBase, parameters, path + simulationName + "_0", mpiWorldSize);
@@ -417,7 +417,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     }
     else
     {
-        Restart::readRestart( dataBase, path + simulationName + "_" + std::to_string( restartIter ) + "_rank_" + std::to_string(rank), startIter );
+        GksGpu::Restart::readRestart( dataBase, path + simulationName + "_" + std::to_string( restartIter ) + "_rank_" + std::to_string(rank), startIter );
 
         if (rank == 0) writeVtkXMLParallelSummaryFile( dataBase, parameters, path + simulationName + "_" + std::to_string( restartIter ) + "_restart", mpiWorldSize );
 
@@ -428,7 +428,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     dataBase->copyDataHostToDevice();
 
-    Initializer::initializeDataUpdate(dataBase);
+    GksGpu::Initializer::initializeDataUpdate(dataBase);
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -436,12 +436,12 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    CupsAnalyzer cupsAnalyzer( dataBase, true, 300.0 );
+    GksGpu::CupsAnalyzer cupsAnalyzer( dataBase, true, 300.0 );
 
-    ConvergenceAnalyzer convergenceAnalyzer( dataBase );
+    GksGpu::ConvergenceAnalyzer convergenceAnalyzer( dataBase );
 
     //auto turbulenceAnalyzer = std::make_shared<TurbulenceAnalyzer>( dataBase, 0 );
-    auto turbulenceAnalyzer = std::make_shared<TurbulenceAnalyzer>( dataBase, 50000 );
+    auto turbulenceAnalyzer = std::make_shared<GksGpu::TurbulenceAnalyzer>( dataBase, 50000 );
 
     turbulenceAnalyzer->collect_UU = true;
     turbulenceAnalyzer->collect_VV = true;
@@ -470,8 +470,8 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
     //    if( subDomainBox->isInside(  0.485, y,  0.3*H ) ) pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(  0.485, y,  0.3*H ), 'W', 10000 );
     //}
 
-    HeatFluxAnalyzer heatFluxAnalyzerPZ(dataBase, bcPZ, 100, 100, lambdaHot, lambdaCold, L);
-    HeatFluxAnalyzer heatFluxAnalyzerMZ(dataBase, bcMZ, 100, 100, lambdaHot, lambdaCold, L);
+    GksGpu::HeatFluxAnalyzer heatFluxAnalyzerPZ(dataBase, bcPZ, 100, 10000, lambdaHot, lambdaCold, L);
+    GksGpu::HeatFluxAnalyzer heatFluxAnalyzerMZ(dataBase, bcMZ, 100, 10000, lambdaHot, lambdaCold, L);
     //HeatFluxAnalyzer heatFluxAnalyzer(dataBase, bcPZ);
 
     //////////////////////////////////////////////////////////////////////////
@@ -480,16 +480,7 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
 
     for( uint iter = startIter + 1; iter <= 100000000; iter++ )
     {
-        TimeStepping::nestedTimeStep(dataBase, parameters, 0);
-
-        if( iter % 10000 == 0 )
-        {
-            dataBase->copyDataDeviceToHost();
-
-            if( rank == 0 ) writeVtkXMLParallelSummaryFile( dataBase, parameters, path + simulationName + "_" + std::to_string( iter ), mpiWorldSize );
-
-            writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
-        }
+        GksGpu::TimeStepping::nestedTimeStep(dataBase, parameters, 0);
 
         cupsAnalyzer.run( iter, parameters.dt );
 
@@ -500,21 +491,33 @@ void simulation( std::string path, std::string simulationName, bool fine, uint r
         if(rankZ == 1) heatFluxAnalyzerPZ.run( iter, parameters );
         if(rankZ == 0) heatFluxAnalyzerMZ.run( iter, parameters );
 
+        if( iter % 10000 == 0 )
+        //if( iter % 25 == 0 )
+        {
+            dataBase->copyDataDeviceToHost();
+
+            if( rank == 0 ) writeVtkXMLParallelSummaryFile( dataBase, parameters, path + simulationName + "_" + std::to_string( iter ), mpiWorldSize );
+
+            writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
+
+            if(rankZ == 1) heatFluxAnalyzerPZ.writeToFile( path + simulationName + "_Nu_top_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
+            if(rankZ == 0) heatFluxAnalyzerMZ.writeToFile( path + simulationName + "_Nu_bot_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
+        }
+
         //pointTimeSeriesCollector->run(iter, parameters);
 
         if( iter > 50000 && iter % 10000 == 0 )
-        //if(iter % 1000 == 0)
         {
             turbulenceAnalyzer->download();
-
+        
             if( rank == 0 ) writeTurbulenceVtkXMLParallelSummaryFile( dataBase, turbulenceAnalyzer, parameters, path + simulationName + "_Turbulence_" + std::to_string( iter ), mpiWorldSize );
-
+        
             writeTurbulenceVtkXML( dataBase, turbulenceAnalyzer, 0, path + simulationName + "_Turbulence_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
         }
 
-        if( iter > 50000 && iter % 10000 == 0 )
+        if( iter % 10000 == 0 )
         {
-            Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ) + "_rank_" + std::to_string(rank), iter );
+            GksGpu::Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ) + "_rank_" + std::to_string(rank), iter );
 
             turbulenceAnalyzer->writeRestartFile( path + simulationName + "_Turbulence_" + std::to_string( iter ) + "_rank_" + std::to_string(rank) );
         }
@@ -549,14 +552,15 @@ int main( int argc, char* argv[])
     int mpiWorldSize = 1;
     MPI_Comm_size(MPI_COMM_WORLD, &mpiWorldSize);
 #else
-    int rank         = MpiUtility::getMpiRankBeforeInit();
-    int mpiWorldSize = MpiUtility::getMpiWorldSizeBeforeInit();
+    int rank         = GksGpu::MpiUtility::getMpiRankBeforeInit();
+    int mpiWorldSize = GksGpu::MpiUtility::getMpiWorldSizeBeforeInit();
 #endif
 
     //////////////////////////////////////////////////////////////////////////
 
 #ifdef _WIN32
-    std::string path( "F:/Work/Computations/out/RayleighBenardMultiGPU/" );
+    std::string path( "F:/Work/Computations/out/RayleighBenardMultiGPU/test/" );
+    //std::string path( "F:/Work/Computations/out/RayleighBenardMultiGPU/" );
 #else
     std::string path( "out/" );
 #endif
@@ -579,7 +583,7 @@ int main( int argc, char* argv[])
     //////////////////////////////////////////////////////////////////////////
 
     // Important: for Cuda-Aware MPI the device must be set before MPI_Init()
-    int deviceCount = CudaUtility::getCudaDeviceCount();
+    int deviceCount = GksGpu::CudaUtility::getCudaDeviceCount();
 
     if(deviceCount == 0)
     {
@@ -588,7 +592,7 @@ int main( int argc, char* argv[])
         *logging::out << logging::Logger::WARNING << msg.str(); msg.str("");
     }
 
-    CudaUtility::setCudaDevice( rank % deviceCount );
+    GksGpu::CudaUtility::setCudaDevice( rank % deviceCount );
 
     //////////////////////////////////////////////////////////////////////////
 
-- 
GitLab