diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 1d2d5677e7014277659fd12e0e2348a5894a7845..b2bf4853b3d822ad308131c69c1de38200d8ae5e 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -73,7 +73,7 @@ Simulation::Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemo
     : para(para), cudaMemoryManager(memoryManager), communicator(communicator), kernelFactory(std::make_unique<KernelFactoryImp>()),
       preProcessorFactory(std::make_shared<PreProcessorFactoryImp>()), dataWriter(std::make_unique<FileWriter>())
 {
-	init(gridProvider, bcFactory);
+    init(gridProvider, bcFactory);
 }
 
 void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFactory)
@@ -104,43 +104,21 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
 
     restart_object = std::make_shared<ASCIIRestartObject>();
     //////////////////////////////////////////////////////////////////////////
-    output.setName(para->getFName() + StringUtil::toString<int>(para->getMyProcessID()) + ".log");
-    if (para->getMyProcessID() == 0)
-        output.setConsoleOut(true);
-    output.clearLogFile();
-    //////////////////////////////////////////////////////////////////////////
     // CUDA streams
     if (para->getUseStreams()) {
         para->getStreamManager()->launchStreams(2u);
         para->getStreamManager()->createCudaEvents();
     }
     //////////////////////////////////////////////////////////////////////////
-    //
-    // output << para->getNeedInterface().at(0) << "\n";
-    // output << para->getNeedInterface().at(1) << "\n";
-    // output << para->getNeedInterface().at(2) << "\n";
-    // output << para->getNeedInterface().at(3) << "\n";
-    // output << para->getNeedInterface().at(4) << "\n";
-    // output << para->getNeedInterface().at(5) << "\n";
-    //////////////////////////////////////////////////////////////////////////
-    // output << "      \t GridX \t GridY \t GridZ \t DistX \t DistY \t DistZ\n";
-    // for (int testout=0; testout<=para->getMaxLevel();testout++)
-    //{
-    //   output << "Level " << testout << ":  " << para->getGridX().at(testout) << " \t " <<
-    //   para->getGridY().at(testout) << " \t " << para->getGridZ().at(testout) << " \t " <<
-    //   para->getDistX().at(testout) << " \t " << para->getDistY().at(testout) << " \t " <<
-    //   para->getDistZ().at(testout) << " \n";
-    //}
-    //////////////////////////////////////////////////////////////////////////
-    output << "LB_Modell:  D3Q" << para->getD3Qxx() << "\n";
-    output << "Re:         " << para->getRe() << "\n";
-    output << "vis_ratio:  " << para->getViscosityRatio() << "\n";
-    output << "u0_ratio:   " << para->getVelocityRatio() << "\n";
-    output << "delta_rho:  " << para->getDensityRatio() << "\n";
-    output << "QuadricLimiters:  " << para->getQuadricLimitersHost()[0] << "\t" << para->getQuadricLimitersHost()[1]
-           << "\t" << para->getQuadricLimitersHost()[2] << "\n";
+    VF_LOG_INFO("LB_Modell:       D3Q{}", para->getD3Qxx());
+    VF_LOG_INFO("Re:              {}", para->getRe());
+    VF_LOG_INFO("vis_ratio:       {}", para->getViscosityRatio());
+    VF_LOG_INFO("u0_ratio:        {}", para->getVelocityRatio());
+    VF_LOG_INFO("delta_rho:       {}", para->getDensityRatio());
+    VF_LOG_INFO("QuadricLimiters: {}, \t{}, \t{}", para->getQuadricLimitersHost()[0],
+                para->getQuadricLimitersHost()[1], para->getQuadricLimitersHost()[2]);
     if (para->getUseAMD())
-        output << "AMD SGS model:  " << para->getSGSConstant() << "\n";
+        VF_LOG_INFO("AMD SGS model:  {}", para->getSGSConstant());
     //////////////////////////////////////////////////////////////////////////
 
     /////////////////////////////////////////////////////////////////////////
@@ -159,20 +137,18 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     //////////////////////////////////////////////////////////////////////////
     // Kernel init
     //////////////////////////////////////////////////////////////////////////
-    output << "make Kernels  "
-           << "\n";
+    VF_LOG_INFO("make Kernels");
     kernels = kernelFactory->makeKernels(para);
 
-    output << "make AD Kernels  "
-           << "\n";
-    if (para->getDiffOn())
+    if (para->getDiffOn()) {
+        VF_LOG_INFO("make AD Kernels");
         adKernels = kernelFactory->makeAdvDifKernels(para);
+    }
 
     //////////////////////////////////////////////////////////////////////////
     // PreProcessor init
     //////////////////////////////////////////////////////////////////////////
-    output << "make Preprocessors  "
-           << "\n";
+    VF_LOG_INFO("make Preprocessors");
     std::vector<PreProcessorType> preProTypes = kernels.at(0)->getPreProcessorTypes();
     preProcessor = preProcessorFactory->makePreProcessor(preProTypes, para);
 
@@ -218,8 +194,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     // Median
     //////////////////////////////////////////////////////////////////////////
     if (para->getCalcMedian()) {
-        output << "alloc Calculation for Mean Values  "
-               << "\n";
+        VF_LOG_INFO("alloc Calculation for Mean Values");
         if (para->getDiffOn())
             allocMedianAD(para.get(), cudaMemoryManager.get());
         else
@@ -230,8 +205,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     // Turbulence Intensity
     //////////////////////////////////////////////////////////////////////////
     if (para->getCalcTurbulenceIntensity()) {
-        output << "alloc arrays for calculating Turbulence Intensity  "
-               << "\n";
+        VF_LOG_INFO("alloc arrays for calculating Turbulence Intensity");
         allocTurbulenceIntensity(para.get(), cudaMemoryManager.get());
     }
 
@@ -255,112 +229,103 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     // MeasurePoints
     //////////////////////////////////////////////////////////////////////////
     if (para->getUseMeasurePoints()) {
-        output << "read measure points...";
+        VF_LOG_INFO("read measure points");
         readMeasurePoints(para.get(), cudaMemoryManager.get());
-        output << "done.\n";
     }
 
     //////////////////////////////////////////////////////////////////////////
     // Porous Media
     //////////////////////////////////////////////////////////////////////////
     if (para->getSimulatePorousMedia()) {
-        output << "define area(s) of porous media...";
+        VF_LOG_INFO("define area(s) of porous media");
         porousMedia();
         kernelFactory->setPorousMedia(pm);
-        output << "done.\n";
     }
 
     //////////////////////////////////////////////////////////////////////////
     // enSightGold
     //////////////////////////////////////////////////////////////////////////
     // excludeGridInterfaceNodesForMirror(para, 7);
-    ////output << "print case file...";
+    ////VF_LOG_INFO("print case file...");
     // printCaseFile(para);
-    ////output << "done.\n";
-    ////output << "print geo file...";
+    ////VF_LOG_INFO("print geo file...");
     // printGeoFile(para, true);  //true for binary
-    ////output << "done.\n";
+    ////VF_LOG_INFO("done.");
 
     //////////////////////////////////////////////////////////////////////////
     // Forcing
     //////////////////////////////////////////////////////////////////////////
     ////allocVeloForForcing(para);
-    // output << "new object forceCalulator  " << "\n";
+    // VF_LOG_INFO("new object forceCalulator");
     // forceCalculator = std::make_shared<ForceCalculations>(para.get());
 
     //////////////////////////////////////////////////////////////////////////
-    // output << "define the Grid..." ;
+    // VF_LOG_INFO("define the Grid...");
     // defineGrid(para, communicator);
     ////allocateMemory();
-    // output << "done.\n";
+    // VF_LOG_INFO("done.");
 
-    output << "init lattice...";
+    VF_LOG_INFO("init lattice...");
     initLattice(para, preProcessor, cudaMemoryManager);
-    output << "done.\n";
+    VF_LOG_INFO("done");
 
-    // output << "set geo for Q...\n" ;
+    // VF_LOG_INFO("set geo for Q...\n");
     // setGeoForQ();
-    // output << "done.\n";
 
     // if (maxlevel>1)
     //{
-    // output << "find Qs...\n" ;
+    // VF_LOG_INFO("find Qs...");
     // findQ27(para);
-    // output << "done.\n";
+    // VF_LOG_INFO("done.");
     //}
 
     // if (para->getDiffOn()==true)
     //{
-    //   output << "define TempBC...\n" ;
+    //   VF_LOG_INFO("define TempBC...");
     //   findTempSim(para);
-    //   output << "done.\n";
 
-    //   output << "define TempVelBC...\n" ;
+    //   VF_LOG_INFO("define TempVelBC...");
     //   findTempVelSim(para);
-    //   output << "done.\n";
 
-    //   output << "define TempPressBC...\n" ;
+    //   VF_LOG_INFO("define TempPressBC...");
     //   findTempPressSim(para);
-    //   output << "done.\n";
+    //   VF_LOG_INFO("done.");
     //}
 
-    // output << "find Qs-BC...\n" ;
+    // VF_LOG_INFO("find Qs-BC...");
     // findBC27(para);
-    // output << "done.\n";
 
-    // output << "find Press-BC...\n" ;
+    // VF_LOG_INFO("find Press-BC...");
     // findPressQShip(para);
-    // output << "done.\n";
+    // VF_LOG_INFO("done.");
 
     //////////////////////////////////////////////////////////////////////////
     // find indices of corner nodes for multiGPU communication
     //////////////////////////////////////////////////////////////////////////
     if (para->getDevices().size() > 2) {
-        output << "Find indices of edge nodes for multiGPU communication ...";
+        VF_LOG_INFO("Find indices of edge nodes for multiGPU communication");
         vf::gpu::findEdgeNodesCommMultiGPU(*para);
-        output << "done.\n";
     }
     //////////////////////////////////////////////////////////////////////////
     // Memory alloc for CheckPoint / Restart
     //////////////////////////////////////////////////////////////////////////
     if (para->getDoCheckPoint() || para->getDoRestart()) {
-        output << "Alloc Memory for CheckPoint / Restart...";
+        VF_LOG_INFO("Alloc Memory for CheckPoint / Restart");
         for (int lev = para->getCoarse(); lev <= para->getFine(); lev++) {
             cudaMemoryManager->cudaAllocFsForCheckPointAndRestart(lev);
         }
-        output << "done.\n";
     }
 
     //////////////////////////////////////////////////////////////////////////
     // Restart
     //////////////////////////////////////////////////////////////////////////
     if (para->getDoRestart()) {
-        output << "Restart...\n...get the Object...\n";
+        VF_LOG_INFO("Restart...\n...get the Object...");
 
         const auto name = getFileName(para->getFName(), para->getTimeDoRestart(), para->getMyProcessID());
         restart_object->deserialize(name, para);
 
-        output << "...copy Memory for Restart...\n";
+        VF_LOG_INFO("...copy Memory for Restart...");
         for (int lev = para->getCoarse(); lev <= para->getFine(); lev++) {
             //////////////////////////////////////////////////////////////////////////
             cudaMemoryManager->cudaCopyFsForRestart(lev);
@@ -377,7 +342,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
             // test...should not work...and does not
             // para->getEvenOrOdd(lev)==false;
         }
-        output << "done.\n";
+        VF_LOG_INFO("done.");
     }
 
     //////////////////////////////////////////////////////////////////////////
@@ -467,27 +432,27 @@ void Simulation::run()
    para->setAngularVelocity(AngularVelocity);
    for (int i = 0; i<= para->getMaxLevel(); i++)
    {
-	   para->getParD(i)->deltaPhi = (real)(para->getAngularVelocity()/(pow(2.,i)));
+       para->getParD(i)->deltaPhi = (real)(para->getAngularVelocity()/(pow(2.,i)));
    }
    //////////////////////////////////////////////////////////////////////////
 
    t_prev = para->getTimeCalcMedStart();
 
-	Timer* averageTimer = new Timer("Average performance");
-	averageTimer->startTimer();
+    Timer* averageTimer = new Timer("Average performance");
+    averageTimer->startTimer();
 
-	////////////////////////////////////////////////////////////////////////////////
-	// Time loop
-	////////////////////////////////////////////////////////////////////////////////
-	for(timestep=para->getTimestepStart();timestep<=para->getTimestepEnd();timestep++)
-	{
+    ////////////////////////////////////////////////////////////////////////////////
+    // Time loop
+    ////////////////////////////////////////////////////////////////////////////////
+    for(timestep=para->getTimestepStart();timestep<=para->getTimestepEnd();timestep++)
+    {
         this->updateGrid27->updateGrid(0, timestep);
 
-	    ////////////////////////////////////////////////////////////////////////////////
-	    //Particles
-	    ////////////////////////////////////////////////////////////////////////////////
-	    if (para->getCalcParticles()) propagateParticles(para.get(), timestep);
-	    ////////////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////////////
+        //Particles
+        ////////////////////////////////////////////////////////////////////////////////
+        if (para->getCalcParticles()) propagateParticles(para.get(), timestep);
+        ////////////////////////////////////////////////////////////////////////////////
 
 
 
@@ -495,14 +460,14 @@ void Simulation::run()
         ////////////////////////////////////////////////////////////////////////////////
         // run Analyzers for kinetic energy and enstrophy for TGV in 3D
         // these analyzers only work on level 0
-	    ////////////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////////////
         if (this->kineticEnergyAnalyzer || this->enstrophyAnalyzer) {
             updateGrid27->exchangeData(0);
         }
 
-	    if( this->kineticEnergyAnalyzer ) this->kineticEnergyAnalyzer->run(timestep);
-	    if( this->enstrophyAnalyzer     ) this->enstrophyAnalyzer->run(timestep);
-	    ////////////////////////////////////////////////////////////////////////////////
+        if( this->kineticEnergyAnalyzer ) this->kineticEnergyAnalyzer->run(timestep);
+        if( this->enstrophyAnalyzer     ) this->enstrophyAnalyzer->run(timestep);
+        ////////////////////////////////////////////////////////////////////////////////
 
 
 
@@ -514,62 +479,62 @@ void Simulation::run()
         {
           for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
           {
-        	  //CalcMedSP27(para->getParD(lev)->vx_SP_Med,
-        			//	  para->getParD(lev)->vy_SP_Med,
-        			//	  para->getParD(lev)->vz_SP_Med,
-        			//	  para->getParD(lev)->rho_SP_Med,
-        			//	  para->getParD(lev)->press_SP_Med,
-        			//	  para->getParD(lev)->geoSP,
-        			//	  para->getParD(lev)->neighborX_SP,
-        			//	  para->getParD(lev)->neighborY_SP,
-        			//	  para->getParD(lev)->neighborZ_SP,
-        			//	  para->getParD(lev)->size_Mat_SP,
-        			//	  para->getParD(lev)->numberofthreads,
-        			//	  para->getParD(lev)->d0SP.f[0],
-        			//	  para->getParD(lev)->evenOrOdd);
-        	  //getLastCudaError("CalcMacSP27 execution failed");
-
-        	  CalcMedCompSP27(para->getParD(lev)->vx_SP_Med,
-        					  para->getParD(lev)->vy_SP_Med,
-        					  para->getParD(lev)->vz_SP_Med,
-        					  para->getParD(lev)->rho_SP_Med,
-        					  para->getParD(lev)->press_SP_Med,
-        					  para->getParD(lev)->typeOfGridNode,
-        					  para->getParD(lev)->neighborX,
-        					  para->getParD(lev)->neighborY,
-        					  para->getParD(lev)->neighborZ,
-        					  para->getParD(lev)->numberOfNodes,
-        					  para->getParD(lev)->numberofthreads,
-        					  para->getParD(lev)->distributions.f[0],
-        					  para->getParD(lev)->isEvenTimestep);
-        	  getLastCudaError("CalcMacMedCompSP27 execution failed");
+              //CalcMedSP27(para->getParD(lev)->vx_SP_Med,
+                    //      para->getParD(lev)->vy_SP_Med,
+                    //      para->getParD(lev)->vz_SP_Med,
+                    //      para->getParD(lev)->rho_SP_Med,
+                    //      para->getParD(lev)->press_SP_Med,
+                    //      para->getParD(lev)->geoSP,
+                    //      para->getParD(lev)->neighborX_SP,
+                    //      para->getParD(lev)->neighborY_SP,
+                    //      para->getParD(lev)->neighborZ_SP,
+                    //      para->getParD(lev)->size_Mat_SP,
+                    //      para->getParD(lev)->numberofthreads,
+                    //      para->getParD(lev)->d0SP.f[0],
+                    //      para->getParD(lev)->evenOrOdd);
+              //getLastCudaError("CalcMacSP27 execution failed");
+
+              CalcMedCompSP27(para->getParD(lev)->vx_SP_Med,
+                              para->getParD(lev)->vy_SP_Med,
+                              para->getParD(lev)->vz_SP_Med,
+                              para->getParD(lev)->rho_SP_Med,
+                              para->getParD(lev)->press_SP_Med,
+                              para->getParD(lev)->typeOfGridNode,
+                              para->getParD(lev)->neighborX,
+                              para->getParD(lev)->neighborY,
+                              para->getParD(lev)->neighborZ,
+                              para->getParD(lev)->numberOfNodes,
+                              para->getParD(lev)->numberofthreads,
+                              para->getParD(lev)->distributions.f[0],
+                              para->getParD(lev)->isEvenTimestep);
+              getLastCudaError("CalcMacMedCompSP27 execution failed");
 
           }
         }
 
-		if (para->getCalcTurbulenceIntensity()) {
+        if (para->getCalcTurbulenceIntensity()) {
             for (int lev = para->getCoarse(); lev <= para->getFine(); lev++) {
-				CalcTurbulenceIntensityDevice(
-				    para->getParD(lev)->vxx,
-				    para->getParD(lev)->vyy,
-				    para->getParD(lev)->vzz,
-				    para->getParD(lev)->vxy,
-				    para->getParD(lev)->vxz,
-				    para->getParD(lev)->vyz,
-				    para->getParD(lev)->vx_mean,
-				    para->getParD(lev)->vy_mean,
-				    para->getParD(lev)->vz_mean,
-				    para->getParD(lev)->distributions.f[0],
-				    para->getParD(lev)->typeOfGridNode,
-				    para->getParD(lev)->neighborX,
-				    para->getParD(lev)->neighborY,
-				    para->getParD(lev)->neighborZ,
-				    para->getParD(lev)->numberOfNodes,
-				    para->getParD(lev)->isEvenTimestep,
-				    para->getParD(lev)->numberofthreads
-				);
-			}
-		}
+                CalcTurbulenceIntensityDevice(
+                    para->getParD(lev)->vxx,
+                    para->getParD(lev)->vyy,
+                    para->getParD(lev)->vzz,
+                    para->getParD(lev)->vxy,
+                    para->getParD(lev)->vxz,
+                    para->getParD(lev)->vyz,
+                    para->getParD(lev)->vx_mean,
+                    para->getParD(lev)->vy_mean,
+                    para->getParD(lev)->vz_mean,
+                    para->getParD(lev)->distributions.f[0],
+                    para->getParD(lev)->typeOfGridNode,
+                    para->getParD(lev)->neighborX,
+                    para->getParD(lev)->neighborY,
+                    para->getParD(lev)->neighborZ,
+                    para->getParD(lev)->numberOfNodes,
+                    para->getParD(lev)->isEvenTimestep,
+                    para->getParD(lev)->numberofthreads
+                );
+            }
+        }
         ////////////////////////////////////////////////////////////////////////////////
 
 
@@ -580,27 +545,27 @@ void Simulation::run()
         ////////////////////////////////////////////////////////////////////////////////
         if(para->getDoCheckPoint() && para->getTimeDoCheckPoint()>0 && timestep%para->getTimeDoCheckPoint()==0 && timestep>0 && !para->overWritingRestart(timestep))
         {
-			averageTimer->stopTimer();
+            averageTimer->stopTimer();
             //////////////////////////////////////////////////////////////////////////
 
             if( para->getDoCheckPoint() )
             {
-                output << "Copy data for CheckPoint t=" << timestep << "...\n";
+                VF_LOG_INFO("Copy data for CheckPoint t = {}....", timestep);
 
                 for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
                 {
                     cudaMemoryManager->cudaCopyFsForCheckPoint(lev);
                 }
 
-                output << "Write data for CheckPoint t=" << timestep << "...";
+                VF_LOG_INFO("Write data for CheckPoint t = {}...", timestep);
 
-				const auto name = getFileName(para->getFName(), timestep, para->getMyProcessID());
-				restart_object->serialize(name, para);
+                const auto name = getFileName(para->getFName(), timestep, para->getMyProcessID());
+                restart_object->serialize(name, para);
 
-                output << "\n done\n";
+                VF_LOG_INFO("done");
             }
             //////////////////////////////////////////////////////////////////////////
-			averageTimer->startTimer();
+            averageTimer->startTimer();
         }
         //////////////////////////////////////////////////////////////////////////////
 
@@ -619,13 +584,13 @@ void Simulation::run()
                 unsigned int valuesPerClockCycle = (unsigned int)(para->getclockCycleForMP() / para->getTimestepForMP());
                 for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
                 {
-                    //output << "start level = " << lev << "\n";
-                    LBCalcMeasurePoints27(  para->getParD(lev)->VxMP,			para->getParD(lev)->VyMP,			para->getParD(lev)->VzMP,
-                    				        para->getParD(lev)->RhoMP,		    para->getParD(lev)->kMP,			para->getParD(lev)->numberOfPointskMP,
-                    				        valuesPerClockCycle,				t_MP,								para->getParD(lev)->typeOfGridNode,
-                    				        para->getParD(lev)->neighborX,   para->getParD(lev)->neighborY,	para->getParD(lev)->neighborZ,
-                    				        para->getParD(lev)->numberOfNodes,	para->getParD(lev)->distributions.f[0],		para->getParD(lev)->numberofthreads,
-                    				        para->getParD(lev)->isEvenTimestep);
+                    // VF_LOG_INFO("start level = {}", lev);
+                    LBCalcMeasurePoints27(  para->getParD(lev)->VxMP,            para->getParD(lev)->VyMP,                 para->getParD(lev)->VzMP,
+                                            para->getParD(lev)->RhoMP,           para->getParD(lev)->kMP,                  para->getParD(lev)->numberOfPointskMP,
+                                            valuesPerClockCycle,                 t_MP,                                     para->getParD(lev)->typeOfGridNode,
+                                            para->getParD(lev)->neighborX,       para->getParD(lev)->neighborY,            para->getParD(lev)->neighborZ,
+                                            para->getParD(lev)->numberOfNodes,   para->getParD(lev)->distributions.f[0],   para->getParD(lev)->numberofthreads,
+                                            para->getParD(lev)->isEvenTimestep);
                 }
                 t_MP++;
             }
@@ -637,7 +602,7 @@ void Simulation::run()
                 {
                     cudaMemoryManager->cudaCopyMeasurePointsToHost(lev);
                     para->copyMeasurePointsArrayToVector(lev);
-                    output << "\n Write MeasurePoints at level = " << lev << " and timestep = " << timestep << "\n";
+                    VF_LOG_INFO("Write MeasurePoints at level = {} and timestep = {}", lev, timestep);
                     for (int j = 0; j < (int)para->getParH(lev)->MP.size(); j++)
                     {
                         MeasurePointWriter::writeMeasurePoints(para.get(), lev, j, timestep);
@@ -658,40 +623,40 @@ void Simulation::run()
         if (para->getDiffOn() && para->getCalcPlaneConc())
         {
             PlaneConcThS27( para->getParD(0)->ConcPlaneIn,
-            		       para->getParD(0)->cpTopIndex,
-            		       para->getParD(0)->numberOfPointsCpTop,
-            		       para->getParD(0)->typeOfGridNode,
-            		       para->getParD(0)->neighborX,
-            		       para->getParD(0)->neighborY,
-            		       para->getParD(0)->neighborZ,
-            		       para->getParD(0)->numberOfNodes,
-            		       para->getParD(0)->numberofthreads,
-            		       para->getParD(0)->distributionsAD27.f[0],
-            		       para->getParD(0)->isEvenTimestep);
+                           para->getParD(0)->cpTopIndex,
+                           para->getParD(0)->numberOfPointsCpTop,
+                           para->getParD(0)->typeOfGridNode,
+                           para->getParD(0)->neighborX,
+                           para->getParD(0)->neighborY,
+                           para->getParD(0)->neighborZ,
+                           para->getParD(0)->numberOfNodes,
+                           para->getParD(0)->numberofthreads,
+                           para->getParD(0)->distributionsAD27.f[0],
+                           para->getParD(0)->isEvenTimestep);
             getLastCudaError("PlaneConcThS27 execution failed");
             PlaneConcThS27( para->getParD(0)->ConcPlaneOut1,
-            		        para->getParD(0)->cpBottomIndex,
-            		        para->getParD(0)->numberOfPointsCpBottom,
-            		        para->getParD(0)->typeOfGridNode,
-            		        para->getParD(0)->neighborX,
-            		        para->getParD(0)->neighborY,
-            		        para->getParD(0)->neighborZ,
-            		        para->getParD(0)->numberOfNodes,
-            		        para->getParD(0)->numberofthreads,
-            		        para->getParD(0)->distributionsAD27.f[0],
-            		        para->getParD(0)->isEvenTimestep);
+                            para->getParD(0)->cpBottomIndex,
+                            para->getParD(0)->numberOfPointsCpBottom,
+                            para->getParD(0)->typeOfGridNode,
+                            para->getParD(0)->neighborX,
+                            para->getParD(0)->neighborY,
+                            para->getParD(0)->neighborZ,
+                            para->getParD(0)->numberOfNodes,
+                            para->getParD(0)->numberofthreads,
+                            para->getParD(0)->distributionsAD27.f[0],
+                            para->getParD(0)->isEvenTimestep);
             getLastCudaError("PlaneConcThS27 execution failed");
             PlaneConcThS27( para->getParD(0)->ConcPlaneOut2,
-            		        para->getParD(0)->pressureBC.kN,
-            		        para->getParD(0)->pressureBC.numberOfBCnodes,
-            		        para->getParD(0)->typeOfGridNode,
-            		        para->getParD(0)->neighborX,
-            		        para->getParD(0)->neighborY,
-            		        para->getParD(0)->neighborZ,
-            		        para->getParD(0)->numberOfNodes,
-            		        para->getParD(0)->numberofthreads,
-            		        para->getParD(0)->distributionsAD27.f[0],
-            		        para->getParD(0)->isEvenTimestep);
+                            para->getParD(0)->pressureBC.kN,
+                            para->getParD(0)->pressureBC.numberOfBCnodes,
+                            para->getParD(0)->typeOfGridNode,
+                            para->getParD(0)->neighborX,
+                            para->getParD(0)->neighborY,
+                            para->getParD(0)->neighborZ,
+                            para->getParD(0)->numberOfNodes,
+                            para->getParD(0)->numberofthreads,
+                            para->getParD(0)->distributionsAD27.f[0],
+                            para->getParD(0)->isEvenTimestep);
             getLastCudaError("PlaneConcThS27 execution failed");
             //////////////////////////////////////////////////////////////////////////////////
             ////Calculation of concentration at the plane
@@ -703,29 +668,29 @@ void Simulation::run()
 
 
 
-	  ////////////////////////////////////////////////////////////////////////////////
+      ////////////////////////////////////////////////////////////////////////////////
       // File IO
       ////////////////////////////////////////////////////////////////////////////////
       //communicator->startTimer();
       if(para->getTimestepOut()>0 && timestep%para->getTimestepOut()==0 && timestep>para->getTimestepStartOut())
       {
-		  //////////////////////////////////////////////////////////////////////////////////
-		  //if (para->getParD(0)->evenOrOdd==true)  para->getParD(0)->evenOrOdd=false;
-		  //else                                    para->getParD(0)->evenOrOdd=true;
-		  //////////////////////////////////////////////////////////////////////////////////
+          //////////////////////////////////////////////////////////////////////////////////
+          //if (para->getParD(0)->evenOrOdd==true)  para->getParD(0)->evenOrOdd=false;
+          //else                                    para->getParD(0)->evenOrOdd=true;
+          //////////////////////////////////////////////////////////////////////////////////
 
-		//////////////////////////////////////////////////////////////////////////
-		averageTimer->stopTimer();
-		averageTimer->outputPerformance(timestep, para.get(), communicator);
-		//////////////////////////////////////////////////////////////////////////
+        //////////////////////////////////////////////////////////////////////////
+        averageTimer->stopTimer();
+        averageTimer->outputPerformance(timestep, para.get(), communicator);
+        //////////////////////////////////////////////////////////////////////////
 
          if( para->getPrintFiles() )
          {
             VF_LOG_INFO("Write files t = {} ...", timestep);
             for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
             {
-		        //////////////////////////////////////////////////////////////////////////
-		        //exchange data for valid post process
+                //////////////////////////////////////////////////////////////////////////
+                //exchange data for valid post process
                 updateGrid27->exchangeData(lev);
                 //////////////////////////////////////////////////////////////////////////
                //if (para->getD3Qxx()==19)
@@ -736,26 +701,26 @@ void Simulation::run()
                //}
                //else if (para->getD3Qxx()==27)
                //{
-				   //if (para->getCalcMedian() && ((int)t > para->getTimeCalcMedStart()) && ((int)t <= para->getTimeCalcMedEnd()))
-				   //{
-					  // unsigned int tdiff = t - t_prev;
-					  // CalcMacMedSP27(para->getParD(lev)->vx_SP_Med,
-				   //					  para->getParD(lev)->vy_SP_Med,
-				   //					  para->getParD(lev)->vz_SP_Med,
-				   //					  para->getParD(lev)->rho_SP_Med,
-				   //					  para->getParD(lev)->press_SP_Med,
-				   //					  para->getParD(lev)->geoSP,
-				   //					  para->getParD(lev)->neighborX_SP,
-				   //					  para->getParD(lev)->neighborY_SP,
-				   //					  para->getParD(lev)->neighborZ_SP,
-				   //					  tdiff,
-				   //					  para->getParD(lev)->size_Mat_SP,
-				   //					  para->getParD(lev)->numberofthreads,
-				   //					  para->getParD(lev)->evenOrOdd);
-					  // getLastCudaError("CalcMacMedSP27 execution failed");
-				   //}
-
-				   //CalcMacSP27(para->getParD(lev)->vx_SP,
+                   //if (para->getCalcMedian() && ((int)t > para->getTimeCalcMedStart()) && ((int)t <= para->getTimeCalcMedEnd()))
+                   //{
+                      // unsigned int tdiff = t - t_prev;
+                      // CalcMacMedSP27(para->getParD(lev)->vx_SP_Med,
+                   //                      para->getParD(lev)->vy_SP_Med,
+                   //                      para->getParD(lev)->vz_SP_Med,
+                   //                      para->getParD(lev)->rho_SP_Med,
+                   //                      para->getParD(lev)->press_SP_Med,
+                   //                      para->getParD(lev)->geoSP,
+                   //                      para->getParD(lev)->neighborX_SP,
+                   //                      para->getParD(lev)->neighborY_SP,
+                   //                      para->getParD(lev)->neighborZ_SP,
+                   //                      tdiff,
+                   //                      para->getParD(lev)->size_Mat_SP,
+                   //                      para->getParD(lev)->numberofthreads,
+                   //                      para->getParD(lev)->evenOrOdd);
+                      // getLastCudaError("CalcMacMedSP27 execution failed");
+                   //}
+
+                   //CalcMacSP27(para->getParD(lev)->vx_SP,
        //                        para->getParD(lev)->vy_SP,
        //                        para->getParD(lev)->vz_SP,
        //                        para->getParD(lev)->rho,
@@ -771,75 +736,75 @@ void Simulation::run()
        //            getLastCudaError("CalcMacSP27 execution failed");
 
 
-				   CalcMacCompSP27(para->getParD(lev)->velocityX,
-								   para->getParD(lev)->velocityY,
-								   para->getParD(lev)->velocityZ,
-								   para->getParD(lev)->rho,
-								   para->getParD(lev)->pressure,
-								   para->getParD(lev)->typeOfGridNode,
-								   para->getParD(lev)->neighborX,
-								   para->getParD(lev)->neighborY,
-								   para->getParD(lev)->neighborZ,
-								   para->getParD(lev)->numberOfNodes,
-								   para->getParD(lev)->numberofthreads,
-								   para->getParD(lev)->distributions.f[0],
-								   para->getParD(lev)->isEvenTimestep);
+                   CalcMacCompSP27(para->getParD(lev)->velocityX,
+                                   para->getParD(lev)->velocityY,
+                                   para->getParD(lev)->velocityZ,
+                                   para->getParD(lev)->rho,
+                                   para->getParD(lev)->pressure,
+                                   para->getParD(lev)->typeOfGridNode,
+                                   para->getParD(lev)->neighborX,
+                                   para->getParD(lev)->neighborY,
+                                   para->getParD(lev)->neighborZ,
+                                   para->getParD(lev)->numberOfNodes,
+                                   para->getParD(lev)->numberofthreads,
+                                   para->getParD(lev)->distributions.f[0],
+                                   para->getParD(lev)->isEvenTimestep);
                    getLastCudaError("CalcMacSP27 execution failed");
 
-				// // overwrite with wall nodes
-				//    SetOutputWallVelocitySP27(  para->getParD(lev)->numberofthreads,
-				// 							   para->getParD(lev)->velocityX,
-				// 							   para->getParD(lev)->velocityY,
-				// 							   para->getParD(lev)->velocityZ,
-				// 							   para->getParD(lev)->geometryBC.Vx,
-				// 							   para->getParD(lev)->geometryBC.Vy,
-				// 							   para->getParD(lev)->geometryBC.Vz,
-				// 							   para->getParD(lev)->geometryBC.numberOfBCnodes,
-				// 							   para->getParD(lev)->geometryBC.k,
-				// 							   para->getParD(lev)->rho,
-				// 							   para->getParD(lev)->pressure,
-				// 							   para->getParD(lev)->typeOfGridNode,
-				// 							   para->getParD(lev)->neighborX,
-				// 							   para->getParD(lev)->neighborY,
-				// 							   para->getParD(lev)->neighborZ,
-				// 							   para->getParD(lev)->size_Mat,
-				// 							   para->getParD(lev)->distributions.f[0],
-				// 							   para->getParD(lev)->isEvenTimestep);
+                // // overwrite with wall nodes
+                //    SetOutputWallVelocitySP27(  para->getParD(lev)->numberofthreads,
+                //                                para->getParD(lev)->velocityX,
+                //                                para->getParD(lev)->velocityY,
+                //                                para->getParD(lev)->velocityZ,
+                //                                para->getParD(lev)->geometryBC.Vx,
+                //                                para->getParD(lev)->geometryBC.Vy,
+                //                                para->getParD(lev)->geometryBC.Vz,
+                //                                para->getParD(lev)->geometryBC.numberOfBCnodes,
+                //                                para->getParD(lev)->geometryBC.k,
+                //                                para->getParD(lev)->rho,
+                //                                para->getParD(lev)->pressure,
+                //                                para->getParD(lev)->typeOfGridNode,
+                //                                para->getParD(lev)->neighborX,
+                //                                para->getParD(lev)->neighborY,
+                //                                para->getParD(lev)->neighborZ,
+                //                                para->getParD(lev)->size_Mat,
+                //                                para->getParD(lev)->distributions.f[0],
+                //                                para->getParD(lev)->isEvenTimestep);
                 //   getLastCudaError("SetOutputWallVelocitySP27 execution failed");
 
-   				//    SetOutputWallVelocitySP27(  para->getParD(lev)->numberofthreads,
-				// 							   para->getParD(lev)->velocityX,
-				// 							   para->getParD(lev)->velocityY,
-				// 							   para->getParD(lev)->velocityZ,
-				// 							   para->getParD(lev)->velocityBC.Vx,
-				// 							   para->getParD(lev)->velocityBC.Vy,
-				// 							   para->getParD(lev)->velocityBC.Vz,
-				// 							   para->getParD(lev)->velocityBC.numberOfBCnodes,
-				// 							   para->getParD(lev)->velocityBC.k,
-				// 							   para->getParD(lev)->rho,
-				// 							   para->getParD(lev)->pressure,
-				// 							   para->getParD(lev)->typeOfGridNode,
-				// 							   para->getParD(lev)->neighborX,
-				// 							   para->getParD(lev)->neighborY,
-				// 							   para->getParD(lev)->neighborZ,
-				// 							   para->getParD(lev)->size_Mat,
-				// 							   para->getParD(lev)->distributions.f[0],
-				// 							   para->getParD(lev)->isEvenTimestep);
+                //    SetOutputWallVelocitySP27(  para->getParD(lev)->numberofthreads,
+                //                                para->getParD(lev)->velocityX,
+                //                                para->getParD(lev)->velocityY,
+                //                                para->getParD(lev)->velocityZ,
+                //                                para->getParD(lev)->velocityBC.Vx,
+                //                                para->getParD(lev)->velocityBC.Vy,
+                //                                para->getParD(lev)->velocityBC.Vz,
+                //                                para->getParD(lev)->velocityBC.numberOfBCnodes,
+                //                                para->getParD(lev)->velocityBC.k,
+                //                                para->getParD(lev)->rho,
+                //                                para->getParD(lev)->pressure,
+                //                                para->getParD(lev)->typeOfGridNode,
+                //                                para->getParD(lev)->neighborX,
+                //                                para->getParD(lev)->neighborY,
+                //                                para->getParD(lev)->neighborZ,
+                //                                para->getParD(lev)->size_Mat,
+                //                                para->getParD(lev)->distributions.f[0],
+                //                                para->getParD(lev)->isEvenTimestep);
                 //   getLastCudaError("SetOutputWallVelocitySP27 execution failed");
 
-				 //}
+                 //}
 
-				   cudaMemoryManager->cudaCopyPrint(lev);
-			   if (para->getCalcMedian())
-			   {
-				   cudaMemoryManager->cudaCopyMedianPrint(lev);
-			   }
+                   cudaMemoryManager->cudaCopyPrint(lev);
+               if (para->getCalcMedian())
+               {
+                   cudaMemoryManager->cudaCopyMedianPrint(lev);
+               }
 
-			   //////////////////////////////////////////////////////////////////////////
+               //////////////////////////////////////////////////////////////////////////
                //TODO: implement flag to write ASCII data
-			   if (para->getWriteVeloASCIIfiles())
-				   VeloASCIIWriter::writeVelocitiesAsTXT(para.get(), lev, timestep);
-			   //////////////////////////////////////////////////////////////////////////
+               if (para->getWriteVeloASCIIfiles())
+                   VeloASCIIWriter::writeVelocitiesAsTXT(para.get(), lev, timestep);
+               //////////////////////////////////////////////////////////////////////////
                if( this->kineticEnergyAnalyzer || this->enstrophyAnalyzer )
                {
                    std::string fname = para->getFName() + "_ID_" + StringUtil::toString<int>(para->getMyProcessID()) + "_t_" + StringUtil::toString<int>(timestep);
@@ -847,10 +812,10 @@ void Simulation::run()
                    if (this->kineticEnergyAnalyzer) this->kineticEnergyAnalyzer->writeToFile(fname);
                    if (this->enstrophyAnalyzer)     this->enstrophyAnalyzer->writeToFile(fname);
                }
-			   //////////////////////////////////////////////////////////////////////////
+               //////////////////////////////////////////////////////////////////////////
 
 
-			   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+               ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
                if (para->getDiffOn()==true)
                {
                   if (para->getDiffMod() == 7)
@@ -880,459 +845,459 @@ void Simulation::run()
                                     para->getParD(lev)->isEvenTimestep);
                   }
 
-				  cudaMemoryManager->cudaCopyConcentrationDeviceToHost(lev);
+                  cudaMemoryManager->cudaCopyConcentrationDeviceToHost(lev);
                   //cudaMemoryCopy(para->getParH(lev)->Conc, para->getParD(lev)->Conc,  para->getParH(lev)->mem_size_real_SP , cudaMemcpyDeviceToHost);
                }
                ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			   ////print cp
-			   //if ((para->getParH(lev)->cpTop.size() > 0) && (t > para->getTStartOut()))
-			   //{
-				  // printCpTopIntermediateStep(para, t, lev);
-			   //}
-			   ////////////////////////////////////////////////////////////////////////////////
-			   //MeasurePointWriter::writeSpacialAverageForXZSlices(para, lev, t);
-			   ////////////////////////////////////////////////////////////////////////////////
-			   //MeasurePointWriter::writeTestAcousticXY(para, lev, t);
-			   //MeasurePointWriter::writeTestAcousticYZ(para, lev, t);
-			   //MeasurePointWriter::writeTestAcousticXZ(para, lev, t);
-			   ////////////////////////////////////////////////////////////////////////
-			}
-
-			//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			////test print press mirror
-			//if (t > para->getTStartOut())
-			//{
-			//	////////////////////////////////////////////////////////////////////////////////
-			//	//Level 7
-			//	CalcCPtop27(para->getParD(7)->d0SP.f[0],
-			//		para->getParD(7)->cpTopIndex,
-			//		para->getParD(7)->numberOfPointsCpTop,
-			//		para->getParD(7)->cpPressTop,
-			//		para->getParD(7)->neighborX_SP,
-			//		para->getParD(7)->neighborY_SP,
-			//		para->getParD(7)->neighborZ_SP,
-			//		para->getParD(7)->size_Mat_SP,
-			//		para->getParD(7)->evenOrOdd,
-			//		para->getParD(7)->numberofthreads);
-			//	//////////////////////////////////////////////////////////////////////////////////
-			//	calcPressForMirror(para, 7);
-			//	////////////////////////////////////////////////////////////////////////////////
-			//	//Level 8
-			//	CalcCPtop27(para->getParD(8)->d0SP.f[0],
-			//		para->getParD(8)->cpTopIndex,
-			//		para->getParD(8)->numberOfPointsCpTop,
-			//		para->getParD(8)->cpPressTop,
-			//		para->getParD(8)->neighborX_SP,
-			//		para->getParD(8)->neighborY_SP,
-			//		para->getParD(8)->neighborZ_SP,
-			//		para->getParD(8)->size_Mat_SP,
-			//		para->getParD(8)->evenOrOdd,
-			//		para->getParD(8)->numberofthreads);
-			//	//////////////////////////////////////////////////////////////////////////////////
-			//	calcPressForMirror(para, 8);
-			//	////////////////////////////////////////////////////////////////////////////////
-			//	//print press mirror
-			//	printScalars(para, false);
-			//	////////////////////////////////////////////////////////////////////////////////
-			//}
-			//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-			//t_prev = t;
-
-			//////////////////////////////////////////////////////////////////////////
-			////Data Analysis
-			////AnalysisData::writeAnalysisData(para, t);
-			//AnalysisData::writeAnalysisDataX(para, t);
-			//AnalysisData::writeAnalysisDataZ(para, t);
-			//////////////////////////////////////////////////////////////////////////
+               ////print cp
+               //if ((para->getParH(lev)->cpTop.size() > 0) && (t > para->getTStartOut()))
+               //{
+                  // printCpTopIntermediateStep(para, t, lev);
+               //}
+               ////////////////////////////////////////////////////////////////////////////////
+               //MeasurePointWriter::writeSpacialAverageForXZSlices(para, lev, t);
+               ////////////////////////////////////////////////////////////////////////////////
+               //MeasurePointWriter::writeTestAcousticXY(para, lev, t);
+               //MeasurePointWriter::writeTestAcousticYZ(para, lev, t);
+               //MeasurePointWriter::writeTestAcousticXZ(para, lev, t);
+               ////////////////////////////////////////////////////////////////////////
+            }
+
+            //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+            ////test print press mirror
+            //if (t > para->getTStartOut())
+            //{
+            //    ////////////////////////////////////////////////////////////////////////////////
+            //    //Level 7
+            //    CalcCPtop27(para->getParD(7)->d0SP.f[0],
+            //        para->getParD(7)->cpTopIndex,
+            //        para->getParD(7)->numberOfPointsCpTop,
+            //        para->getParD(7)->cpPressTop,
+            //        para->getParD(7)->neighborX_SP,
+            //        para->getParD(7)->neighborY_SP,
+            //        para->getParD(7)->neighborZ_SP,
+            //        para->getParD(7)->size_Mat_SP,
+            //        para->getParD(7)->evenOrOdd,
+            //        para->getParD(7)->numberofthreads);
+            //    //////////////////////////////////////////////////////////////////////////////////
+            //    calcPressForMirror(para, 7);
+            //    ////////////////////////////////////////////////////////////////////////////////
+            //    //Level 8
+            //    CalcCPtop27(para->getParD(8)->d0SP.f[0],
+            //        para->getParD(8)->cpTopIndex,
+            //        para->getParD(8)->numberOfPointsCpTop,
+            //        para->getParD(8)->cpPressTop,
+            //        para->getParD(8)->neighborX_SP,
+            //        para->getParD(8)->neighborY_SP,
+            //        para->getParD(8)->neighborZ_SP,
+            //        para->getParD(8)->size_Mat_SP,
+            //        para->getParD(8)->evenOrOdd,
+            //        para->getParD(8)->numberofthreads);
+            //    //////////////////////////////////////////////////////////////////////////////////
+            //    calcPressForMirror(para, 8);
+            //    ////////////////////////////////////////////////////////////////////////////////
+            //    //print press mirror
+            //    printScalars(para, false);
+            //    ////////////////////////////////////////////////////////////////////////////////
+            //}
+            //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+            //t_prev = t;
+
+            //////////////////////////////////////////////////////////////////////////
+            ////Data Analysis
+            ////AnalysisData::writeAnalysisData(para, t);
+            //AnalysisData::writeAnalysisDataX(para, t);
+            //AnalysisData::writeAnalysisDataZ(para, t);
+            //////////////////////////////////////////////////////////////////////////
 
             ////////////////////////////////////////////////////////////////////////
             //pressure difference
             ////////////////////////////////////////////////////////////////////////
-			   //if (para->getMyID() == para->getPressInID())       calcPressure(para,  "in", 0);
-			   //else if (para->getMyID() == para->getPressOutID()) calcPressure(para, "out", 0);
+               //if (para->getMyID() == para->getPressInID())       calcPressure(para,  "in", 0);
+               //else if (para->getMyID() == para->getPressOutID()) calcPressure(para, "out", 0);
             ////////////////////////////////////////////////////////////////////////
             //flow rate
             ////////////////////////////////////////////////////////////////////////
-		      //calcFlowRate(para, 0);
+              //calcFlowRate(para, 0);
+            ////////////////////////////////////////////////////////////////////////
+
+            ////////////////////////////////////////////////////////////////////////
+            //calculate 2nd, 3rd and higher order moments
+            ////////////////////////////////////////////////////////////////////////
+            if (para->getCalc2ndOrderMoments())  calc2ndMoments(para.get(), cudaMemoryManager.get());
+            if (para->getCalc3rdOrderMoments())  calc3rdMoments(para.get(), cudaMemoryManager.get());
+            if (para->getCalcHighOrderMoments()) calcHigherOrderMoments(para.get(), cudaMemoryManager.get());
             ////////////////////////////////////////////////////////////////////////
 
-			////////////////////////////////////////////////////////////////////////
-			//calculate 2nd, 3rd and higher order moments
-			////////////////////////////////////////////////////////////////////////
-			if (para->getCalc2ndOrderMoments())  calc2ndMoments(para.get(), cudaMemoryManager.get());
-			if (para->getCalc3rdOrderMoments())  calc3rdMoments(para.get(), cudaMemoryManager.get());
-			if (para->getCalcHighOrderMoments()) calcHigherOrderMoments(para.get(), cudaMemoryManager.get());
-			////////////////////////////////////////////////////////////////////////
-
-			////////////////////////////////////////////////////////////////////////
-			//calculate median on host
-			////////////////////////////////////////////////////////////////////////
-			if (para->getCalcMedian() && ((int)timestep > para->getTimeCalcMedStart()) && ((int)timestep <= para->getTimeCalcMedEnd()) && ((timestep%(unsigned int)para->getclockCycleForMP())==0))
-			{
-				unsigned int tdiff = timestep - t_prev;
-				calcMedian(para.get(), tdiff);
-
-				/////////////////////////////////
-				//added for incremental averaging
-				t_prev = timestep;
-				resetMedian(para.get());
-				/////////////////////////////////
-			}
+            ////////////////////////////////////////////////////////////////////////
+            //calculate median on host
+            ////////////////////////////////////////////////////////////////////////
+            if (para->getCalcMedian() && ((int)timestep > para->getTimeCalcMedStart()) && ((int)timestep <= para->getTimeCalcMedEnd()) && ((timestep%(unsigned int)para->getclockCycleForMP())==0))
+            {
+                unsigned int tdiff = timestep - t_prev;
+                calcMedian(para.get(), tdiff);
+
+                /////////////////////////////////
+                //added for incremental averaging
+                t_prev = timestep;
+                resetMedian(para.get());
+                /////////////////////////////////
+            }
             if (para->getCalcTurbulenceIntensity())
-			{
+            {
                 uint t_diff = timestep - t_turbulenceIntensity;
                 calcTurbulenceIntensity(para.get(), cudaMemoryManager.get(), t_diff);
                 //writeAllTiDatafToFile(para.get(), t);
             }
-			////////////////////////////////////////////////////////////////////////
-			dataWriter->writeTimestep(para, timestep);
-			////////////////////////////////////////////////////////////////////////
+            ////////////////////////////////////////////////////////////////////////
+            dataWriter->writeTimestep(para, timestep);
+            ////////////////////////////////////////////////////////////////////////
             if (para->getCalcTurbulenceIntensity()) {
                 t_turbulenceIntensity = timestep;
                 resetVelocityFluctuationsAndMeans(para.get(), cudaMemoryManager.get());
             }
-			////////////////////////////////////////////////////////////////////////
+            ////////////////////////////////////////////////////////////////////////
             if (para->getCalcDragLift()) printDragLift(para.get(), cudaMemoryManager.get(), timestep);
-			////////////////////////////////////////////////////////////////////////
-			if (para->getCalcParticles()) copyAndPrintParticles(para.get(), cudaMemoryManager.get(), timestep, false);
-			////////////////////////////////////////////////////////////////////////
+            ////////////////////////////////////////////////////////////////////////
+            if (para->getCalcParticles()) copyAndPrintParticles(para.get(), cudaMemoryManager.get(), timestep, false);
+            ////////////////////////////////////////////////////////////////////////
             VF_LOG_INFO("... done");
-			////////////////////////////////////////////////////////////////////////
+            ////////////////////////////////////////////////////////////////////////
          }
 
-		////////////////////////////////////////////////////////////////////////
-		averageTimer->startTimer();
+        ////////////////////////////////////////////////////////////////////////
+        averageTimer->startTimer();
       }
-	}
-
-	/////////////////////////////////////////////////////////////////////////
-
-	////////////////////////////////////////////////////////////////////////////////
-	//printDragLift(para);
-	////////////////////////////////////////////////////////////////////////////////
-
-	////////////////////////////////////////////////////////////////////////////////
-	if (para->getDiffOn()==true) printPlaneConc(para.get(), cudaMemoryManager.get());
-	////////////////////////////////////////////////////////////////////////////////
-
-	////////////////////////////////////////////////////////////////////////////////
-	////for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
-	////{
-	////	if (para->getParH(lev)->cpTop.size() > 0)
-	////	{
-	////		printCpTop(para, lev);
-	////	}
-	////}
-	//for (int lev = 7; lev <= 8; lev++)
-	//{
-	//	printCpTop(para, lev);
-	//}
-	////printCpTop(para);
-	////printCpBottom(para);
-	////printCpBottom2(para);
-	////////////////////////////////////////////////////////////////////////////////
+    }
+
+    /////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////
+    //printDragLift(para);
+    ////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////
+    if (para->getDiffOn()==true) printPlaneConc(para.get(), cudaMemoryManager.get());
+    ////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////
+    ////for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
+    ////{
+    ////    if (para->getParH(lev)->cpTop.size() > 0)
+    ////    {
+    ////        printCpTop(para, lev);
+    ////    }
+    ////}
+    //for (int lev = 7; lev <= 8; lev++)
+    //{
+    //    printCpTop(para, lev);
+    //}
+    ////printCpTop(para);
+    ////printCpBottom(para);
+    ////printCpBottom2(para);
+    ////////////////////////////////////////////////////////////////////////////////
 
  //  //////////////////////////////////////////////////////////////////////////
  //  //Copy Measure Values
-	//for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
-	//{
-	//	output << "\n Copy MeasurePoints at level = " << lev <<"\n";
-	//	para->cudaCopyMeasurePointsToHost(lev);
-	//	para->copyMeasurePointsArrayToVector(lev);
-	//	output << "\n Write MeasurePoints at level = " << lev <<"\n";
-	//	for(int j = 0; j < (int)para->getParH(lev)->MP.size(); j++)
-	//	{
-	//		MeasurePointWriter::writeMeasurePoints(para, lev, j, 0);
-	//	}
-	//}
+    //for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
+    //{
+    //    VF_LOG_INFO("Copy MeasurePoints at level = {}", lev);
+    //    para->cudaCopyMeasurePointsToHost(lev);
+    //    para->copyMeasurePointsArrayToVector(lev);
+    //    VF_LOG_INFO("Write MeasurePoints at level = {}", lev);
+    //    for(int j = 0; j < (int)para->getParH(lev)->MP.size(); j++)
+    //    {
+    //        MeasurePointWriter::writeMeasurePoints(para, lev, j, 0);
+    //    }
+    //}
  //  //////////////////////////////////////////////////////////////////////////
 }
 
 void Simulation::porousMedia()
 {
-	double porosity, darcySI, forchheimerSI;
-	double dxLBM = 0.00390625;
-	double dtLBM = 0.00000658;
-	unsigned int level, geo;
-	double startX, startY, startZ, endX, endY, endZ;
-	//////////////////////////////////////////////////////////////////////////
-
-	////////////////////////////////////////////////////////////////////////////
-	////Test = porous media 0
-	//porosity = 0.7;
-	//darcySI = 137.36; //[1/s]
-	//forchheimerSI = 1037.8; //[1/m]
-	//level = para->getFine();
-	//geo = GEO_PM_0;
-	//startX = 20.0;
-	//startY =  0.0;
-	//startZ =  0.0;
-	//endX = 40.0;
-	//endY = 22.0;
-	//endZ = 22.0;
-	//pm[0] = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
-	//pm[0]->setStartCoordinates(startX, startY, startZ);
-	//pm[0]->setEndCoordinates(endX, endY, endZ);
-	//pm[0]->setResistanceLBM();
-	//definePMarea(pm[0]);
-	////////////////////////////////////////////////////////////////////////////
-
-	//////////////////////////////////////////////////////////////////////////
-	//Kondensator = porous media 0
-	porosity = 0.7;
-	darcySI = 137.36; //[1/s]
-	forchheimerSI = 1037.8; //[1/m]
-	level = para->getFine();
-	geo = GEO_PM_0;
-	startX = -0.715882;
-	startY = -0.260942;
-	startZ = -0.031321;
-	endX = -0.692484;
-	endY =  0.277833;
-	endZ =  0.360379;
-	pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
-	int n = (int)pm.size() - 1;
-	pm.at(n)->setStartCoordinates(startX, startY, startZ);
-	pm.at(n)->setEndCoordinates(endX, endY, endZ);
-	pm.at(n)->setResistanceLBM();
-	definePMarea(pm.at(n));
-	//////////////////////////////////////////////////////////////////////////
-
-	//////////////////////////////////////////////////////////////////////////
-	//NT-Kuehler = porous media 1
-	porosity = 0.6;
-	darcySI = 149.98; //[1/s]
-	forchheimerSI = 960.57; //[1/m]
-	level = para->getFine();
-	geo = GEO_PM_1;
-	startX = -0.696146;
-	startY = -0.32426;
-	startZ = -0.0421345;
-	endX = -0.651847;
-	endY =  0.324822;
-	endZ =  0.057098;
-	pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
-	n = (int)pm.size() - 1;
-	pm.at(n)->setStartCoordinates(startX, startY, startZ);
-	pm.at(n)->setEndCoordinates(endX, endY, endZ);
-	pm.at(n)->setResistanceLBM();
-	definePMarea(pm.at(n));
-	//////////////////////////////////////////////////////////////////////////
-
-	//////////////////////////////////////////////////////////////////////////
-	//Wasserkuehler = porous media 2
-	porosity = 0.6;
-	darcySI = 148.69; //[1/s]
-	forchheimerSI = 629.45; //[1/m]
-	level = para->getFine();
-	geo = GEO_PM_2;
-	startX = -0.692681;
-	startY = -0.324954;
-	startZ = 0.0789429;
-	endX = -0.657262;
-	endY =  0.32538;
-	endZ =  0.400974;
-	pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
-	n = (int)pm.size() - 1;
-	pm.at(n)->setStartCoordinates(startX, startY, startZ);
-	pm.at(n)->setEndCoordinates(endX, endY, endZ);
-	pm.at(n)->setResistanceLBM();
-	definePMarea(pm.at(n));
-	//////////////////////////////////////////////////////////////////////////
+    double porosity, darcySI, forchheimerSI;
+    double dxLBM = 0.00390625;
+    double dtLBM = 0.00000658;
+    unsigned int level, geo;
+    double startX, startY, startZ, endX, endY, endZ;
+    //////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////
+    ////Test = porous media 0
+    //porosity = 0.7;
+    //darcySI = 137.36; //[1/s]
+    //forchheimerSI = 1037.8; //[1/m]
+    //level = para->getFine();
+    //geo = GEO_PM_0;
+    //startX = 20.0;
+    //startY =  0.0;
+    //startZ =  0.0;
+    //endX = 40.0;
+    //endY = 22.0;
+    //endZ = 22.0;
+    //pm[0] = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+    //pm[0]->setStartCoordinates(startX, startY, startZ);
+    //pm[0]->setEndCoordinates(endX, endY, endZ);
+    //pm[0]->setResistanceLBM();
+    //definePMarea(pm[0]);
+    ////////////////////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////////////////////
+    //Kondensator = porous media 0
+    porosity = 0.7;
+    darcySI = 137.36; //[1/s]
+    forchheimerSI = 1037.8; //[1/m]
+    level = para->getFine();
+    geo = GEO_PM_0;
+    startX = -0.715882;
+    startY = -0.260942;
+    startZ = -0.031321;
+    endX = -0.692484;
+    endY =  0.277833;
+    endZ =  0.360379;
+    pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
+    int n = (int)pm.size() - 1;
+    pm.at(n)->setStartCoordinates(startX, startY, startZ);
+    pm.at(n)->setEndCoordinates(endX, endY, endZ);
+    pm.at(n)->setResistanceLBM();
+    definePMarea(pm.at(n));
+    //////////////////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////////////////////
+    //NT-Kuehler = porous media 1
+    porosity = 0.6;
+    darcySI = 149.98; //[1/s]
+    forchheimerSI = 960.57; //[1/m]
+    level = para->getFine();
+    geo = GEO_PM_1;
+    startX = -0.696146;
+    startY = -0.32426;
+    startZ = -0.0421345;
+    endX = -0.651847;
+    endY =  0.324822;
+    endZ =  0.057098;
+    pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
+    n = (int)pm.size() - 1;
+    pm.at(n)->setStartCoordinates(startX, startY, startZ);
+    pm.at(n)->setEndCoordinates(endX, endY, endZ);
+    pm.at(n)->setResistanceLBM();
+    definePMarea(pm.at(n));
+    //////////////////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////////////////////
+    //Wasserkuehler = porous media 2
+    porosity = 0.6;
+    darcySI = 148.69; //[1/s]
+    forchheimerSI = 629.45; //[1/m]
+    level = para->getFine();
+    geo = GEO_PM_2;
+    startX = -0.692681;
+    startY = -0.324954;
+    startZ = 0.0789429;
+    endX = -0.657262;
+    endY =  0.32538;
+    endZ =  0.400974;
+    pm.push_back(std::shared_ptr<PorousMedia>(new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level)));
+    n = (int)pm.size() - 1;
+    pm.at(n)->setStartCoordinates(startX, startY, startZ);
+    pm.at(n)->setEndCoordinates(endX, endY, endZ);
+    pm.at(n)->setResistanceLBM();
+    definePMarea(pm.at(n));
+    //////////////////////////////////////////////////////////////////////////
 
 }
 
 void Simulation::definePMarea(std::shared_ptr<PorousMedia>& pMedia)
 {
-	unsigned int counter = 0;
-	unsigned int level = pMedia->getLevelPM();
-	std::vector< unsigned int > nodeIDsPorousMedia;
-	output << "definePMarea....find nodes \n";
-
-	for (unsigned int i = 0; i < para->getParH(level)->numberOfNodes; i++)
-	{
-		if (((para->getParH(level)->coordinateX[i] >= pMedia->getStartX()) && (para->getParH(level)->coordinateX[i] <= pMedia->getEndX())) &&
-			((para->getParH(level)->coordinateY[i] >= pMedia->getStartY()) && (para->getParH(level)->coordinateY[i] <= pMedia->getEndY())) &&
-			((para->getParH(level)->coordinateZ[i] >= pMedia->getStartZ()) && (para->getParH(level)->coordinateZ[i] <= pMedia->getEndZ())) )
-		{
-			if (para->getParH(level)->typeOfGridNode[i] >= GEO_FLUID)
-			{
-				para->getParH(level)->typeOfGridNode[i] = pMedia->getGeoID();
-				nodeIDsPorousMedia.push_back(i);
-				counter++;
-			}
-		}
-	}
-
-	output << "definePMarea....cuda copy SP \n";
-	cudaMemoryManager->cudaCopySP(level);
-	pMedia->setSizePM(counter);
-	output << "definePMarea....cuda alloc PM \n";
-	cudaMemoryManager->cudaAllocPorousMedia(pMedia.get(), level);
-	unsigned int *tpmArrayIDs = pMedia->getHostNodeIDsPM();
-
-	output << "definePMarea....copy vector to array \n";
-	for (unsigned int j = 0; j < pMedia->getSizePM(); j++)
-	{
-		tpmArrayIDs[j] = nodeIDsPorousMedia[j];
-	}
-
-	pMedia->setHostNodeIDsPM(tpmArrayIDs);
-	output << "definePMarea....cuda copy PM \n";
-	cudaMemoryManager->cudaCopyPorousMedia(pMedia.get(), level);
+    unsigned int counter = 0;
+    unsigned int level = pMedia->getLevelPM();
+    std::vector< unsigned int > nodeIDsPorousMedia;
+    VF_LOG_INFO("definePMarea....find nodes");
+
+    for (unsigned int i = 0; i < para->getParH(level)->numberOfNodes; i++)
+    {
+        if (((para->getParH(level)->coordinateX[i] >= pMedia->getStartX()) && (para->getParH(level)->coordinateX[i] <= pMedia->getEndX())) &&
+            ((para->getParH(level)->coordinateY[i] >= pMedia->getStartY()) && (para->getParH(level)->coordinateY[i] <= pMedia->getEndY())) &&
+            ((para->getParH(level)->coordinateZ[i] >= pMedia->getStartZ()) && (para->getParH(level)->coordinateZ[i] <= pMedia->getEndZ())) )
+        {
+            if (para->getParH(level)->typeOfGridNode[i] >= GEO_FLUID)
+            {
+                para->getParH(level)->typeOfGridNode[i] = pMedia->getGeoID();
+                nodeIDsPorousMedia.push_back(i);
+                counter++;
+            }
+        }
+    }
+
+    VF_LOG_INFO("definePMarea....cuda copy SP");
+    cudaMemoryManager->cudaCopySP(level);
+    pMedia->setSizePM(counter);
+    VF_LOG_INFO("definePMarea....cuda alloc PM");
+    cudaMemoryManager->cudaAllocPorousMedia(pMedia.get(), level);
+    unsigned int *tpmArrayIDs = pMedia->getHostNodeIDsPM();
+
+    VF_LOG_INFO("definePMarea....copy vector to array");
+    for (unsigned int j = 0; j < pMedia->getSizePM(); j++)
+    {
+        tpmArrayIDs[j] = nodeIDsPorousMedia[j];
+    }
+
+    pMedia->setHostNodeIDsPM(tpmArrayIDs);
+    VF_LOG_INFO("definePMarea....cuda copy PM");
+    cudaMemoryManager->cudaCopyPorousMedia(pMedia.get(), level);
 }
 
 Simulation::~Simulation()
 {
-	// Cuda Streams
+    // Cuda Streams
     if (para->getUseStreams()) {
         para->getStreamManager()->destroyCudaEvents();
         para->getStreamManager()->terminateStreams();
-	}
+    }
 
-	//CudaFreeHostMemory
+    //CudaFreeHostMemory
     for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
-	{
-		//para->cudaFreeFull(lev);
-		cudaMemoryManager->cudaFreeCoord(lev);
-		cudaMemoryManager->cudaFreeSP(lev);
-		if (para->getCalcMedian())
-		{
-			cudaMemoryManager->cudaFreeMedianSP(lev);
-		}
-		//para->cudaFreeVeloBC(lev);
-		//para->cudaFreeWallBC(lev);
-		//para->cudaFreeVeloBC(lev);
-		//para->cudaFreeInlet(lev);
-		//para->cudaFreeOutlet(lev);
-		//para->cudaFreeGeomBC(lev);
-		//para->cudaFreePress(lev);
-	}
-	if (para->getMaxLevel()>1)
-	{
-		for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
-		{
-			cudaMemoryManager->cudaFreeInterfaceCF(lev);
-			cudaMemoryManager->cudaFreeInterfaceFC(lev);
-			cudaMemoryManager->cudaFreeInterfaceOffCF(lev);
-			cudaMemoryManager->cudaFreeInterfaceOffFC(lev);
-			//para->cudaFreePressX1(lev);
-		}
-	}
-	//para->cudaFreeVeloBC(0); //level = 0
-	//para->cudaFreePressBC();
-	//para->cudaFreeVeloPropeller(para->getFine());
-	//para->cudaFreePressX0(para->getCoarse());
-
-	//////////////////////////////////////////////////////////////////////////
-	//Temp
-	if (para->getDiffOn() == true)
-	{
-		for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
-		{
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->Conc_Full));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->Conc));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->Temp.temp));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->Temp.k));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.temp));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.velo));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.k));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.temp));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.velo));
-			checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.k));
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-
-
-	//////////////////////////////////////////////////////////////////////////
-	//free second order moments
-	if (para->getCalc2ndOrderMoments())
-	{
-		for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
-		{
-			cudaMemoryManager->cudaFree2ndMoments(lev);
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-	//free third order moments
-	if (para->getCalc3rdOrderMoments())
-	{
-		for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
-		{
-			cudaMemoryManager->cudaFree3rdMoments(lev);
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-	//free higher order moments
-	if (para->getCalcHighOrderMoments())
-	{
-		for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
-		{
-			cudaMemoryManager->cudaFreeHigherMoments(lev);
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-
-
-	//////////////////////////////////////////////////////////////////////////
-	//Multi GPU
-	//////////////////////////////////////////////////////////////////////////
-	////1D domain decomposition
-	//if (para->getNumprocs() > 1)
-	//{
-	// for (int lev=para->getCoarse(); lev < para->getFine(); lev++)
-	// {
-	//  for (unsigned int i=0; i < para->getNumberOfProcessNeighbors(lev, "send"); i++)
-	//  {
-	//   para->cudaFreeProcessNeighbor(lev, i);
-	//  }
-	// }
-	//}
-	//////////////////////////////////////////////////////////////////////////
-	//3D domain decomposition
-	if (para->getNumprocs() > 1)
-	{
-		for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
-		{
-			//////////////////////////////////////////////////////////////////////////
-			for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsX(lev, "send"); i++)
-			{
-				cudaMemoryManager->cudaFreeProcessNeighborX(lev, i);
-			}
-			//////////////////////////////////////////////////////////////////////////
-			for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsY(lev, "send"); i++)
-			{
-				cudaMemoryManager->cudaFreeProcessNeighborY(lev, i);
-			}
-			//////////////////////////////////////////////////////////////////////////
-			for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsZ(lev, "send"); i++)
-			{
-				cudaMemoryManager->cudaFreeProcessNeighborZ(lev, i);
-			}
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-	//Normals
-	if (para->getIsGeoNormal()) {
-		for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
-		{
-			cudaMemoryManager->cudaFreeGeomNormals(lev);
-		}
-	}
-	//////////////////////////////////////////////////////////////////////////
-	// Turbulence Intensity
-	if (para->getCalcTurbulenceIntensity()) {
+    {
+        //para->cudaFreeFull(lev);
+        cudaMemoryManager->cudaFreeCoord(lev);
+        cudaMemoryManager->cudaFreeSP(lev);
+        if (para->getCalcMedian())
+        {
+            cudaMemoryManager->cudaFreeMedianSP(lev);
+        }
+        //para->cudaFreeVeloBC(lev);
+        //para->cudaFreeWallBC(lev);
+        //para->cudaFreeVeloBC(lev);
+        //para->cudaFreeInlet(lev);
+        //para->cudaFreeOutlet(lev);
+        //para->cudaFreeGeomBC(lev);
+        //para->cudaFreePress(lev);
+    }
+    if (para->getMaxLevel()>1)
+    {
+        for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
+        {
+            cudaMemoryManager->cudaFreeInterfaceCF(lev);
+            cudaMemoryManager->cudaFreeInterfaceFC(lev);
+            cudaMemoryManager->cudaFreeInterfaceOffCF(lev);
+            cudaMemoryManager->cudaFreeInterfaceOffFC(lev);
+            //para->cudaFreePressX1(lev);
+        }
+    }
+    //para->cudaFreeVeloBC(0); //level = 0
+    //para->cudaFreePressBC();
+    //para->cudaFreeVeloPropeller(para->getFine());
+    //para->cudaFreePressX0(para->getCoarse());
+
+    //////////////////////////////////////////////////////////////////////////
+    //Temp
+    if (para->getDiffOn() == true)
+    {
+        for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
+        {
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->Conc_Full));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->Conc));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->Temp.temp));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->Temp.k));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.temp));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.velo));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempVel.k));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.temp));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.velo));
+            checkCudaErrors(cudaFreeHost(para->getParH(lev)->TempPress.k));
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+
+
+    //////////////////////////////////////////////////////////////////////////
+    //free second order moments
+    if (para->getCalc2ndOrderMoments())
+    {
+        for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
+        {
+            cudaMemoryManager->cudaFree2ndMoments(lev);
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+    //free third order moments
+    if (para->getCalc3rdOrderMoments())
+    {
+        for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
+        {
+            cudaMemoryManager->cudaFree3rdMoments(lev);
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+    //free higher order moments
+    if (para->getCalcHighOrderMoments())
+    {
+        for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
+        {
+            cudaMemoryManager->cudaFreeHigherMoments(lev);
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+
+
+    //////////////////////////////////////////////////////////////////////////
+    //Multi GPU
+    //////////////////////////////////////////////////////////////////////////
+    ////1D domain decomposition
+    //if (para->getNumprocs() > 1)
+    //{
+    // for (int lev=para->getCoarse(); lev < para->getFine(); lev++)
+    // {
+    //  for (unsigned int i=0; i < para->getNumberOfProcessNeighbors(lev, "send"); i++)
+    //  {
+    //   para->cudaFreeProcessNeighbor(lev, i);
+    //  }
+    // }
+    //}
+    //////////////////////////////////////////////////////////////////////////
+    //3D domain decomposition
+    if (para->getNumprocs() > 1)
+    {
+        for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
+        {
+            //////////////////////////////////////////////////////////////////////////
+            for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsX(lev, "send"); i++)
+            {
+                cudaMemoryManager->cudaFreeProcessNeighborX(lev, i);
+            }
+            //////////////////////////////////////////////////////////////////////////
+            for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsY(lev, "send"); i++)
+            {
+                cudaMemoryManager->cudaFreeProcessNeighborY(lev, i);
+            }
+            //////////////////////////////////////////////////////////////////////////
+            for (unsigned int i = 0; i < para->getNumberOfProcessNeighborsZ(lev, "send"); i++)
+            {
+                cudaMemoryManager->cudaFreeProcessNeighborZ(lev, i);
+            }
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+    //Normals
+    if (para->getIsGeoNormal()) {
+        for (int lev = para->getCoarse(); lev < para->getFine(); lev++)
+        {
+            cudaMemoryManager->cudaFreeGeomNormals(lev);
+        }
+    }
+    //////////////////////////////////////////////////////////////////////////
+    // Turbulence Intensity
+    if (para->getCalcTurbulenceIntensity()) {
         cudaFreeTurbulenceIntensityArrays(para.get(), cudaMemoryManager.get());
-	//PreCollisionInteractors
-	for( SPtr<PreCollisionInteractor> actuator: para->getActuators()){
-		actuator->free(para.get(), cudaMemoryManager.get());
-	}
-
-	for( SPtr<PreCollisionInteractor> probe: para->getProbes()){
-		probe->free(para.get(), cudaMemoryManager.get());
-	}
-	//////////////////////////////////////////////////////////////////////////
+    //PreCollisionInteractors
+    for( SPtr<PreCollisionInteractor> actuator: para->getActuators()){
+        actuator->free(para.get(), cudaMemoryManager.get());
+    }
+
+    for( SPtr<PreCollisionInteractor> probe: para->getProbes()){
+        probe->free(para.get(), cudaMemoryManager.get());
+    }
+    //////////////////////////////////////////////////////////////////////////
     }
-}
\ No newline at end of file
+}
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.h b/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
index 98dca98747d0eb6603812f75428242e05edd23e4..3c3954fe909428a0e20f2755cb8befdbf64b1263 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.h
@@ -6,7 +6,6 @@
 
 #include <PointerDefinitions.h>
 
-#include "Output/LogWriter.hpp"
 #include "Utilities/Buffer2D.hpp"
 #include "LBM/LB.h"
 
@@ -66,8 +65,6 @@ private:
 	Buffer2D <int> geo_rbuf_b;
 
 
-	LogWriter output;
-
 	vf::gpu::Communicator& communicator;
     SPtr<Parameter> para;
     std::shared_ptr<DataWriter> dataWriter;
diff --git a/src/gpu/VirtualFluids_GPU/Output/LogWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/LogWriter.hpp
deleted file mode 100644
index cbce9a48cdc8ca36127ab61363fdb9a3dff9dec8..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Output/LogWriter.hpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#ifndef LOGWRITER_H
-#define LOGWRITER_H
-
-#include <iostream>
-#include <fstream>
-//#include <string>
-
-//#include "Utilities/StringUtil.hpp"
-
-
-////////////////////////////////////////////////////////////////////////////////
-class LogWriter
-{
-public:
-   LogWriter()
-   {  
-      consoleOut = false;
-   }
-   LogWriter(std::string fname)
-   {
-      consoleOut = false;
-      this->fname = fname;
-   }
-   void setName(std::string name)
-   {
-      this->fname = name;
-   }
-   void setConsoleOut(bool flag)
-   {
-      consoleOut = flag;
-   }
-   void clearLogFile()
-   {
-      ostr.open(fname.c_str(), std::ios_base::out);
-      if (ostr.bad())
-      {
-         std::string exceptionText = "Error: Output file/directory not found! LogWriter::operator << \n";
-         throw exceptionText;
-      }
-      ostr << "";
-      ostr.close();
-   }
-   template <typename T>
-   LogWriter&  operator << (const T& arg)
-   {
-      ostr.open(fname.c_str(), std::ios_base::app);
-      if (ostr.bad())
-      {
-         //std::cout << "Error: Output file/directory not found! LogWriter::operator <<" << std::endl;
-         //return *this;
-         std::string exceptionText = "Error: Output file/directory not found! LogWriter::operator << \n";
-         throw exceptionText;
-      }
-      ostr << arg;
-      ostr.close();
-      if(consoleOut) std::cout << arg << std::flush;
-      return *this;
-   }
-protected:
-private:
-   std::string fname;
-   std::ofstream ostr;
-   bool consoleOut;
-};
-
-#endif	
diff --git a/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp
index 648c1fb5d5480c12f703f7dfbbd717ec1c515da2..83f0a677b0012153cf079b466a333acc58bda6be 100644
--- a/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp
+++ b/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp
@@ -58,4 +58,4 @@ inline void writeNeighborLinkLinesDebug(Parameter *para)
 
 } // namespace NeighborDebugWriter
 
-#endif
\ No newline at end of file
+#endif
diff --git a/src/lbm/constants/D3Q27.h b/src/lbm/constants/D3Q27.h
index 4cd9a5be8f070e7eaa18c3114075b4d728f3a73a..6a198e926477eff4534108686793400de3a7e042 100644
--- a/src/lbm/constants/D3Q27.h
+++ b/src/lbm/constants/D3Q27.h
@@ -11,35 +11,35 @@ static constexpr int STARTDIR = 0;
 static constexpr int ENDDIR   = 26;
 
 // used in the CPU and the GPU version
-static constexpr int DIR_000 = 0;	 // REST
-static constexpr int DIR_P00 = 1;	 // E
-static constexpr int DIR_M00 = 2;	 // W
-static constexpr int DIR_0P0 = 3;	 // N
-static constexpr int DIR_0M0 = 4;	 // S
-static constexpr int DIR_00P = 5;	 // T
-static constexpr int DIR_00M = 6;	 // B
-
-static constexpr int DIR_PP0 = 7;	 // NE
-static constexpr int DIR_MM0 = 8;	 // SW
-static constexpr int DIR_PM0 = 9;	 // SE
-static constexpr int DIR_MP0 = 10;	 // NW
-static constexpr int DIR_P0P = 11;	 // TE
-static constexpr int DIR_M0M = 12;	 // BW
-static constexpr int DIR_P0M = 13;	 // BE
-static constexpr int DIR_M0P = 14;	 // TW
-static constexpr int DIR_0PP = 15;	 // TN
-static constexpr int DIR_0MM = 16;	 // BS
-static constexpr int DIR_0PM = 17;	 // BN
-static constexpr int DIR_0MP = 18;	 // TS
-
-static constexpr int DIR_PPP = 19;	 // TNE
-static constexpr int DIR_MPP = 20;	 // TNW
-static constexpr int DIR_PMP = 21;	 // TSE
-static constexpr int DIR_MMP = 22;	 // TSW
-static constexpr int DIR_PPM = 23;	 // BNE
-static constexpr int DIR_MPM = 24;	 // BNW
-static constexpr int DIR_PMM = 25;	 // BSE
-static constexpr int DIR_MMM = 26;	 // BSW
+static constexpr int DIR_000 = 0;    // REST
+static constexpr int DIR_P00 = 1;    // E
+static constexpr int DIR_M00 = 2;    // W
+static constexpr int DIR_0P0 = 3;    // N
+static constexpr int DIR_0M0 = 4;    // S
+static constexpr int DIR_00P = 5;    // T
+static constexpr int DIR_00M = 6;    // B
+
+static constexpr int DIR_PP0 = 7;    // NE
+static constexpr int DIR_MM0 = 8;    // SW
+static constexpr int DIR_PM0 = 9;    // SE
+static constexpr int DIR_MP0 = 10;   // NW
+static constexpr int DIR_P0P = 11;   // TE
+static constexpr int DIR_M0M = 12;   // BW
+static constexpr int DIR_P0M = 13;   // BE
+static constexpr int DIR_M0P = 14;   // TW
+static constexpr int DIR_0PP = 15;   // TN
+static constexpr int DIR_0MM = 16;   // BS
+static constexpr int DIR_0PM = 17;   // BN
+static constexpr int DIR_0MP = 18;   // TS
+
+static constexpr int DIR_PPP = 19;   // TNE
+static constexpr int DIR_MPP = 20;   // TNW
+static constexpr int DIR_PMP = 21;   // TSE
+static constexpr int DIR_MMP = 22;   // TSW
+static constexpr int DIR_PPM = 23;   // BNE
+static constexpr int DIR_MPM = 24;   // BNW
+static constexpr int DIR_PMM = 25;   // BSE
+static constexpr int DIR_MMM = 26;   // BSW
 
 struct countersForPointerChasing{
     uint counterInverse;