diff --git a/CMake/cmake_config_files/BOMBADIL.config.cmake b/CMake/cmake_config_files/BOMBADIL.config.cmake
index 9a205596d22cfac18baef3131871b9f73e1b4f9c..c0c5cf2f08b2cc925be248441c99a336160fd1bd 100644
--- a/CMake/cmake_config_files/BOMBADIL.config.cmake
+++ b/CMake/cmake_config_files/BOMBADIL.config.cmake
@@ -29,15 +29,15 @@ set(VTK_DIR "d:/Tools/VTK/build/VTK-8.0.0")
 #################################################################################
 #  METIS  
 #################################################################################
-IF(${USE_METIS})
-  SET(METIS_INCLUDEDIR "d:/Tools/metis-5.1.0/include")
-  SET(METIS_DEBUG_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Debug/metis.lib") 
-  SET(METIS_RELEASE_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Release/metis.lib") 
+#IF(${USE_METIS})
+#  SET(METIS_INCLUDEDIR "d:/Tools/metis-5.1.0/include")
+#  SET(METIS_DEBUG_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Debug/metis.lib") 
+#  SET(METIS_RELEASE_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Release/metis.lib") 
   
   # SET(METIS_INCLUDEDIR "/mnt/d/Tools/metis-5.1.0/include")
   # SET(METIS_DEBUG_LIBRARY "/mnt/d/Tools/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a") 
   # SET(METIS_RELEASE_LIBRARY "/mnt/d/Tools/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a") 
-ENDIF()
+#ENDIF()
 
 #################################################################################
 #  PE  
diff --git a/CMake/cmake_config_files/LISE.config.cmake b/CMake/cmake_config_files/LISE.config.cmake
index 42f10b7dd54854d8b9e4a50616468a7f9292a331..705f02c62f5beb1eb5544af88d9a248c65e684ba 100644
--- a/CMake/cmake_config_files/LISE.config.cmake
+++ b/CMake/cmake_config_files/LISE.config.cmake
@@ -1,6 +1,3 @@
-LIST(APPEND CAB_ADDTIONAL_COMPILER_FLAGS -D__unix__)
-LIST(APPEND CAB_ADDTIONAL_COMPILER_FLAGS -D__UNIX__)
-
 #################################################################################
 #  METIS  
 #################################################################################
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 4bf538b84b0211ada593913a556f10be0d6ed608..a0693d564438ece97e0e9bebc7262bf189a337d4 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -1,12 +1,12 @@
 # Contributing
 
 If you want to contribute to VirtualFluids, your help is very welcome.
-To contribute use a pull request as follows:
+To contribute use a merge request as follows:
 
-### How to make a clean pull request
+### How to make a clean merge request
 
 - Create a personal fork of VirtualFluids.
-- Clone the fork on your local machine. Your remote repo on gitea is called `origin`.
+- Clone the fork on your local machine. Your remote repo on gitlab is called `origin`.
 - Add the original repository as a remote called `upstream`.
 - If you created your fork a while ago be sure to pull upstream changes into your local repository.
 - Create a new branch to work on! Branch from `develop` or `open_source`.
@@ -16,10 +16,10 @@ To contribute use a pull request as follows:
 - Write or adapt tests as needed.
 - Add or change the documentation as needed.
 - Push your branch to your fork on gitea, the remote `origin`.
-- From your fork open a pull request in the correct branch. Target the project's `develop` or `open_source` branch
+- From your fork open a merge request in the correct branch. Target the project's `develop` or `open_source` branch
 - …
-- If we requests further changes just push them to your branch. The PR will be updated automatically.
-- Once the pull request is approved and merged you can pull the changes from `upstream` to your local repo and delete
+- If we requests further changes just push them to your branch. The MR will be updated automatically.
+- Once the merge request is approved and merged you can pull the changes from `upstream` to your local repo and delete
 your extra branch(es).
 
 And last but not least: Always write your commit messages in the present tense. Your commit message should describe what the commit, when applied, does to the code – not what you did to the code.
diff --git a/README.md b/README.md
index 615e4d758e95b4b556efee9f9487beaf08d90d21..dcd9c19d6ed9a2254cb021904c3cfd982ad237e8 100644
--- a/README.md
+++ b/README.md
@@ -11,8 +11,8 @@ VirtualFluids has been used on a variety of platforms:
  - Cygwin
 ### Software Requirements
  
- - [CMake](https://cmake.org/) (minimum version 3.13)
- - C++ compiler with C++11 support, for example gcc 6.3 or Visual C++ 14.0
+ - [CMake](https://cmake.org/) (minimum version 3.15)
+ - C++ compiler with C++14 support
  - [Paraview](https://www.paraview.org/) (most recent version)
 
 with usage of the GPU:  
@@ -62,9 +62,7 @@ The doxygen generated documentation can be found [here](https://git.irmb.bau.tu-
 
 
 ## Known Issues
-If CMake does not find CUDA_CUT_INCLUDE_DIR use and set the correct CUDA Pathes in gpu.cmake in the base directory in lines 35, 36.
-
-If you notice any problems on your platform, please report an gitea issue. 
+If you notice any problems on your platform, please report an issue. 
 
 
 ## Authors
diff --git a/ansible/playbook_cppcheck.yml b/ansible/playbook_cppcheck.yml
deleted file mode 100644
index 942f842450c91d27f3eea08c21060e152c7c2fed..0000000000000000000000000000000000000000
--- a/ansible/playbook_cppcheck.yml
+++ /dev/null
@@ -1,10 +0,0 @@
-- hosts: gitlab_ci_deploy_cppcheck
-  tasks:
-    - name: Create remote cppcheck dir
-      command: mkdir ~/cppcheck
-      ignore_errors: yes
-
-    - name: Synchronize cppcheck with remote
-      synchronize:
-        src: "../html_report"
-        dest: "~/cppcheck"
diff --git a/ansible/playbook_gcov.yml b/ansible/playbook_gcov.yml
deleted file mode 100644
index 45e53b9cd08090750a3e10e7e1e2c9ffa71e0766..0000000000000000000000000000000000000000
--- a/ansible/playbook_gcov.yml
+++ /dev/null
@@ -1,10 +0,0 @@
-- hosts: gitlab_ci_deploy_gcov
-  tasks:
-    - name: Create remote gcov dir
-      command: mkdir ~/gcov
-      ignore_errors: yes
-
-    - name: Synchronize gcov with remote
-      synchronize:
-        src: "../coverage"
-        dest: "~/gcov"
diff --git a/apps/cpu/Multiphase/Multiphase.cfg b/apps/cpu/Multiphase/Multiphase.cfg
index 0d3e5633ea5b4c718eadd5fba04b89abc38497ed..34ef25bbc20475c16e011191f652ed1a1de3a2f4 100644
--- a/apps/cpu/Multiphase/Multiphase.cfg
+++ b/apps/cpu/Multiphase/Multiphase.cfg
@@ -44,5 +44,5 @@ restartStep = 100000
 cpStart = 100000
 cpStep = 100000
 
-outTime = 1
-endTime = 200000000
\ No newline at end of file
+outTime = 100
+endTime = 2000
\ No newline at end of file
diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index 1973ae061c02118c3a08210ece90f1b6905ef81a..f122e90157d8ecbdc7128a156edba407414ad2d0 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -81,7 +81,7 @@ void run(string configname)
 
         SPtr<LBMKernel> kernel;
 
-        kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
+        kernel = SPtr<LBMKernel>(new MultiphaseScratchCumulantLBMKernel());
 
         kernel->setWithForcing(true);
         kernel->setForcingX1(0.0);
@@ -471,7 +471,7 @@ void run(string configname)
         grid->accept(setConnsVisitor);
 
         SPtr<UbScheduler> visSch(new UbScheduler(outTime));
-        SPtr<WriteMacroscopicQuantitiesCoProcessor> pp(new WriteMacroscopicQuantitiesCoProcessor(
+        SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
             grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), conv, comm));
 
         SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
diff --git a/src/basics/Core/Input/ConfigData/ConfigDataImp.cpp b/src/basics/Core/Input/ConfigData/ConfigDataImp.cpp
index 7e0e093610b6f1decde482abddb9ba331213810d..0c53b32829eb989667b93fc7df156525267e930b 100644
--- a/src/basics/Core/Input/ConfigData/ConfigDataImp.cpp
+++ b/src/basics/Core/Input/ConfigData/ConfigDataImp.cpp
@@ -5,86 +5,6 @@ std::shared_ptr<ConfigDataImp> ConfigDataImp::getNewInstance()
     return std::shared_ptr<ConfigDataImp>(new ConfigDataImp());
 }
 
-ConfigDataImp::ConfigDataImp()
-{
-    this->isViscosity            = false;
-    this->isNumberOfDevices      = false;
-    this->isDevices              = false;
-    this->isOutputPath           = false;
-    this->isPrefix               = false;
-    this->isGridPath             = false;
-    this->isPrintOutputFiles     = false;
-    this->isGeometryValues       = false;
-    this->isCalc2ndOrderMoments  = false;
-    this->isCalc3rdOrderMoments  = false;
-    this->isCalcHighOrderMoments = false;
-    this->isReadGeo              = false;
-    this->isCalcMedian           = false;
-    this->isCalcDragLift         = false;
-    this->isCalcCp               = false;
-    this->isConcFile             = false;
-    this->isUseMeasurePoints     = false;
-    this->isUseWale              = false;
-    this->isSimulatePorousMedia  = false;
-    this->isD3Qxx                = false;
-    this->isTEnd                 = false;
-    this->isTOut                 = false;
-    this->isTStartOut            = false;
-    this->isTimeCalcMedStart     = false;
-    this->isTimeCalcMedEnd       = false;
-    this->isPressInID            = false;
-    this->isPressOutID           = false;
-    this->isPressInZ             = false;
-    this->isPressOutZ            = false;
-    this->isDiffOn               = false;
-    this->isDiffMod              = false;
-    this->isDiffusivity          = false;
-    this->isTemperatureInit      = false;
-    this->isTemperatureBC        = false;
-    this->isVelocity             = false;
-    this->isViscosityRatio       = false;
-    this->isVelocityRatio        = false;
-    this->isDensityRatio         = false;
-    this->isPressRatio           = false;
-    this->isRealX                = false;
-    this->isRealY                = false;
-    this->isFactorPressBC        = false;
-    this->isGeometryFileC        = false;
-    this->isGeometryFileM        = false;
-    this->isGeometryFileF        = false;
-    this->isClockCycleForMP      = false;
-    this->isTimestepForMP        = false;
-    this->isForcingX             = false;
-    this->isForcingY             = false;
-    this->isForcingZ             = false;
-    this->isCalcParticles        = false;
-    this->isParticleBasicLevel   = false;
-    this->isParticleInitLevel    = false;
-    this->isNumberOfParticles    = false;
-    this->isNeighborWSB          = false;
-    this->isStartXHotWall        = false;
-    this->isEndXHotWall          = false;
-    this->isPossNeighborFilesX   = false;
-    this->isPossNeighborFilesY   = false;
-    this->isPossNeighborFilesZ   = false;
-    this->isTimeDoCheckPoint     = false;
-    this->isTimeDoRestart        = false;
-    this->isDoCheckPoint         = false;
-    this->isDoRestart            = false;
-    this->isMaxLevel             = false;
-    this->isGridX                = false;
-    this->isGridY                = false;
-    this->isGridZ                = false;
-    this->isDistX                = false;
-    this->isDistY                = false;
-    this->isDistZ                = false;
-    this->isNeedInterface        = false;
-    this->isMainKernel           = false;
-    this->isMultiKernelOn        = false;
-    this->isMultiKernelLevel     = false;
-    this->isMultiKernelName      = false;
-}
-
 real ConfigDataImp::getViscosity() { return this->viscosity; }
 
 uint ConfigDataImp::getNumberOfDevices() { return this->numberOfDevices; }
@@ -633,19 +553,19 @@ void ConfigDataImp::setEndXHotWall(real endXHotWall)
     this->isEndXHotWall = true;
 }
 
-void ConfigDataImp::setPossNeighborFilesX(std::vector<std::string> possNeighborFilesX)
+void ConfigDataImp::setPossNeighborFilesX(const std::vector<std::string> &possNeighborFilesX)
 {
     this->possNeighborFilesX   = possNeighborFilesX;
     this->isPossNeighborFilesX = true;
 }
 
-void ConfigDataImp::setPossNeighborFilesY(std::vector<std::string> possNeighborFilesY)
+void ConfigDataImp::setPossNeighborFilesY(const std::vector<std::string> &possNeighborFilesY)
 {
     this->possNeighborFilesY   = possNeighborFilesY;
     this->isPossNeighborFilesY = true;
 }
 
-void ConfigDataImp::setPossNeighborFilesZ(std::vector<std::string> possNeighborFilesZ)
+void ConfigDataImp::setPossNeighborFilesZ(const std::vector<std::string> &possNeighborFilesZ)
 {
     this->possNeighborFilesZ   = possNeighborFilesZ;
     this->isPossNeighborFilesZ = true;
@@ -681,49 +601,49 @@ void ConfigDataImp::setMaxLevel(uint maxLevel)
     this->isMaxLevel = true;
 }
 
-void ConfigDataImp::setGridX(std::vector<int> gridX)
+void ConfigDataImp::setGridX(const std::vector<int> &gridX)
 {
     this->gridX   = gridX;
     this->isGridX = true;
 }
 
-void ConfigDataImp::setGridY(std::vector<int> gridY)
+void ConfigDataImp::setGridY(const std::vector<int> &gridY)
 {
     this->gridY   = gridY;
     this->isGridY = true;
 }
 
-void ConfigDataImp::setGridZ(std::vector<int> gridZ)
+void ConfigDataImp::setGridZ(const std::vector<int> &gridZ)
 {
     this->gridZ   = gridZ;
     this->isGridZ = true;
 }
 
-void ConfigDataImp::setDistX(std::vector<int> distX)
+void ConfigDataImp::setDistX(const std::vector<int> &distX)
 {
     this->distX   = distX;
     this->isDistX = true;
 }
 
-void ConfigDataImp::setDistY(std::vector<int> distY)
+void ConfigDataImp::setDistY(const std::vector<int> &distY)
 {
     this->distY   = distY;
     this->isDistY = true;
 }
 
-void ConfigDataImp::setDistZ(std::vector<int> distZ)
+void ConfigDataImp::setDistZ(const std::vector<int> &distZ)
 {
     this->distZ   = distZ;
     this->isDistZ = true;
 }
 
-void ConfigDataImp::setNeedInterface(std::vector<bool> needInterface)
+void ConfigDataImp::setNeedInterface(const std::vector<bool> &needInterface)
 {
     this->needInterface   = needInterface;
     this->isNeedInterface = true;
 }
 
-void ConfigDataImp::setMainKernel(std::string mainKernel)
+void ConfigDataImp::setMainKernel(const std::string &mainKernel)
 {
     this->mainKernel   = mainKernel;
     this->isMainKernel = true;
@@ -735,13 +655,13 @@ void ConfigDataImp::setMultiKernelOn(bool multiKernelOn)
     this->isMultiKernelOn = true;
 }
 
-void ConfigDataImp::setMultiKernelLevel(std::vector<int> multiKernelLevel)
+void ConfigDataImp::setMultiKernelLevel(const std::vector<int> &multiKernelLevel)
 {
     this->multiKernelLevel   = multiKernelLevel;
     this->isMultiKernelLevel = true;
 }
 
-void ConfigDataImp::setMultiKernelName(std::vector<std::string> multiKernelName)
+void ConfigDataImp::setMultiKernelName(const std::vector<std::string> &multiKernelName)
 {
     this->multiKernelName   = multiKernelName;
     this->isMultiKernelName = true;
diff --git a/src/basics/Core/Input/ConfigData/ConfigDataImp.h b/src/basics/Core/Input/ConfigData/ConfigDataImp.h
index 34a2a95b919f9ccd23072d0a0608aa269386ef80..511ac4be166c7e43885209f144f0cdac6821c6f1 100644
--- a/src/basics/Core/Input/ConfigData/ConfigDataImp.h
+++ b/src/basics/Core/Input/ConfigData/ConfigDataImp.h
@@ -161,9 +161,9 @@ public:
     void setNumberOfParticles(int numberOfParticles);
     void setStartXHotWall(real startXHotWall);
     void setEndXHotWall(real endXHotWall);
-    void setPossNeighborFilesX(std::vector<std::string> possNeighborFilesX);
-    void setPossNeighborFilesY(std::vector<std::string> possNeighborFilesY);
-    void setPossNeighborFilesZ(std::vector<std::string> possNeighborFilesZ);
+    void setPossNeighborFilesX(const std::vector<std::string> &possNeighborFilesX);
+    void setPossNeighborFilesY(const std::vector<std::string> &possNeighborFilesY);
+    void setPossNeighborFilesZ(const std::vector<std::string> &possNeighborFilesZ);
     // void setPossNeighborFilesX(std::vector<std::string> possNeighborFilesX);
     // void setPossNeighborFilesY(std::vector<std::string> possNeighborFilesY);
     // void setPossNeighborFilesZ(std::vector<std::string> possNeighborFilesZ);
@@ -172,17 +172,17 @@ public:
     void setDoCheckPoint(bool doCheckPoint);
     void setDoRestart(bool doRestart);
     void setMaxLevel(uint maxLevel);
-    void setGridX(std::vector<int> gridX);
-    void setGridY(std::vector<int> gridY);
-    void setGridZ(std::vector<int> gridZ);
-    void setDistX(std::vector<int> distX);
-    void setDistY(std::vector<int> distY);
-    void setDistZ(std::vector<int> distZ);
-    void setNeedInterface(std::vector<bool> needInterface);
-    void setMainKernel(std::string mainKernel);
+    void setGridX(const std::vector<int> &gridX);
+    void setGridY(const std::vector<int> &gridY);
+    void setGridZ(const std::vector<int> &gridZ);
+    void setDistX(const std::vector<int> &distX);
+    void setDistY(const std::vector<int> &distY);
+    void setDistZ(const std::vector<int> &distZ);
+    void setNeedInterface(const std::vector<bool> &needInterface);
+    void setMainKernel(const std::string &mainKernel);
     void setMultiKernelOn(bool multiKernelOn);
-    void setMultiKernelLevel(std::vector<int> multiKernelLevel);
-    void setMultiKernelName(std::vector<std::string> multiKernelName);
+    void setMultiKernelLevel(const std::vector<int> &multiKernelLevel);
+    void setMultiKernelName(const std::vector<std::string> &multiKernelName);
 
     bool isViscosityInConfigFile() override;
     bool isNumberOfDevicesInConfigFile() override;
@@ -270,83 +270,83 @@ public:
     bool isMultiKernelNameInConfigFile() override;
 
 private:
-    ConfigDataImp();
+    ConfigDataImp() = default;
 
-    real viscosity;
-    uint numberOfDevices;
+    real viscosity { 0. };
+    uint numberOfDevices { 0 };
     std::vector<uint> devices;
     std::string outputPath;
     std::string prefix;
     std::string gridPath;
-    bool printOutputFiles;
-    bool geometryValues;
-    bool calc2ndOrderMoments;
-    bool calc3rdOrderMoments;
-    bool calcHighOrderMoments;
-    bool readGeo;
-    bool calcMedian;
-    bool calcDragLift;
-    bool calcCp;
-    bool writeVeloASCIIfiles;
-    bool calcPlaneConc;
-    bool concFile;
-    bool streetVelocityFile;
-    bool useMeasurePoints;
-    bool useWale;
-    bool useInitNeq;
-    bool simulatePorousMedia;
-    uint d3Qxx;
-    uint tEnd;
-    uint tOut;
-    uint tStartOut;
-    uint timeCalcMedStart;
-    uint timeCalcMedEnd;
-    uint pressInID;
-    uint pressOutID;
-    uint pressInZ;
-    uint pressOutZ;
-    bool diffOn;
-    uint diffMod;
-    real diffusivity;
-    real temperatureInit;
-    real temperatureBC;
-    // real viscosity;
-    real velocity;
-    real viscosityRatio;
-    real velocityRatio;
-    real densityRatio;
-    real pressRatio;
-    real realX;
-    real realY;
-    real factorPressBC;
+    bool printOutputFiles { false };
+    bool geometryValues { false };
+    bool calc2ndOrderMoments { false };
+    bool calc3rdOrderMoments { false };
+    bool calcHighOrderMoments { false };
+    bool readGeo { false };
+    bool calcMedian { false };
+    bool calcDragLift { false };
+    bool calcCp { false };
+    bool writeVeloASCIIfiles { false };
+    bool calcPlaneConc { false };
+    bool concFile { false };
+    bool streetVelocityFile { false };
+    bool useMeasurePoints { false };
+    bool useWale { false };
+    bool useInitNeq { false };
+    bool simulatePorousMedia { false };
+    uint d3Qxx { 0 };
+    uint tEnd { 0 };
+    uint tOut { 0 };
+    uint tStartOut { 0 };
+    uint timeCalcMedStart { 0 };
+    uint timeCalcMedEnd { 0 };
+    uint pressInID { 0 };
+    uint pressOutID { 0 };
+    uint pressInZ { 0 };
+    uint pressOutZ { 0 };
+    bool diffOn { false };
+    uint diffMod { 0 };
+    real diffusivity { 0. };
+    real temperatureInit { 0. };
+    real temperatureBC { 0. };
+    // real viscosity { 0 };
+    real velocity { 0. };
+    real viscosityRatio { 0. };
+    real velocityRatio { 0. };
+    real densityRatio { 0. };
+    real pressRatio { 0. };
+    real realX { 0. };
+    real realY { 0. };
+    real factorPressBC { 0. };
     std::string geometryFileC;
     std::string geometryFileM;
     std::string geometryFileF;
-    uint clockCycleForMP;
-    uint timestepForMP;
-    real forcingX;
-    real forcingY;
-    real forcingZ;
-    real quadricLimiterP;
-    real quadricLimiterM;
-    real quadricLimiterD;
-    bool calcParticles;
-    int particleBasicLevel;
-    int particleInitLevel;
-    int numberOfParticles;
-    real startXHotWall;
-    real endXHotWall;
+    uint clockCycleForMP { 0 };
+    uint timestepForMP { 0 };
+    real forcingX { 0. };
+    real forcingY { 0. };
+    real forcingZ { 0. };
+    real quadricLimiterP { 0. };
+    real quadricLimiterM { 0. };
+    real quadricLimiterD { 0. };
+    bool calcParticles { false };
+    int particleBasicLevel { 0 };
+    int particleInitLevel { 0 };
+    int numberOfParticles { 0 };
+    real startXHotWall { 0. };
+    real endXHotWall { 0. };
     std::vector<std::string> possNeighborFilesX;
     std::vector<std::string> possNeighborFilesY;
     std::vector<std::string> possNeighborFilesZ;
     // std::vector<std::string> possNeighborFilesX;
     // std::vector<std::string> possNeighborFilesY;
     // std::vector<std::string> possNeighborFilesZ;
-    int timeDoCheckPoint;
-    int timeDoRestart;
-    bool doCheckPoint;
-    bool doRestart;
-    int maxLevel;
+    int timeDoCheckPoint { 0 };
+    int timeDoRestart { 0 };
+    bool doCheckPoint{ false };
+    bool doRestart{ false };
+    int maxLevel { 0 };
     std::vector<int> gridX;
     std::vector<int> gridY;
     std::vector<int> gridZ;
@@ -355,93 +355,93 @@ private:
     std::vector<int> distZ;
     std::vector<bool> needInterface;
     std::string mainKernel;
-    bool multiKernelOn;
+    bool multiKernelOn{ false };
     std::vector<int> multiKernelLevel;
     std::vector<std::string> multiKernelName;
 
-    bool isViscosity;
-    bool isNumberOfDevices;
-    bool isDevices;
-    bool isOutputPath;
-    bool isPrefix;
-    bool isGridPath;
-    bool isPrintOutputFiles;
-    bool isGeometryValues;
-    bool isCalc2ndOrderMoments;
-    bool isCalc3rdOrderMoments;
-    bool isCalcHighOrderMoments;
-    bool isReadGeo;
-    bool isCalcMedian;
-    bool isCalcDragLift;
-    bool isCalcCp;
-    bool isWriteVeloASCII;
-    bool isCalcPlaneConc;
-    bool isConcFile;
-    bool isStreetVelocityFile;
-    bool isUseMeasurePoints;
-    bool isUseWale;
-    bool isUseInitNeq;
-    bool isSimulatePorousMedia;
-    bool isD3Qxx;
-    bool isTEnd;
-    bool isTOut;
-    bool isTStartOut;
-    bool isTimeCalcMedStart;
-    bool isTimeCalcMedEnd;
-    bool isPressInID;
-    bool isPressOutID;
-    bool isPressInZ;
-    bool isPressOutZ;
-    bool isDiffOn;
-    bool isDiffMod;
-    bool isDiffusivity;
-    bool isTemperatureInit;
-    bool isTemperatureBC;
-    // bool isViscosity;
-    bool isVelocity;
-    bool isViscosityRatio;
-    bool isVelocityRatio;
-    bool isDensityRatio;
-    bool isPressRatio;
-    bool isRealX;
-    bool isRealY;
-    bool isFactorPressBC;
-    bool isGeometryFileC;
-    bool isGeometryFileM;
-    bool isGeometryFileF;
-    bool isClockCycleForMP;
-    bool isTimestepForMP;
-    bool isForcingX;
-    bool isForcingY;
-    bool isForcingZ;
-    bool isQuadricLimiterP;
-    bool isQuadricLimiterM;
-    bool isQuadricLimiterD;
-    bool isCalcParticles;
-    bool isParticleBasicLevel;
-    bool isParticleInitLevel;
-    bool isNumberOfParticles;
-    bool isNeighborWSB;
-    bool isStartXHotWall;
-    bool isEndXHotWall;
-    bool isPossNeighborFilesX;
-    bool isPossNeighborFilesY;
-    bool isPossNeighborFilesZ;
-    bool isTimeDoCheckPoint;
-    bool isTimeDoRestart;
-    bool isDoCheckPoint;
-    bool isDoRestart;
-    bool isMaxLevel;
-    bool isGridX;
-    bool isGridY;
-    bool isGridZ;
-    bool isDistX;
-    bool isDistY;
-    bool isDistZ;
-    bool isNeedInterface;
-    bool isMainKernel;
-    bool isMultiKernelOn;
-    bool isMultiKernelLevel;
-    bool isMultiKernelName;
+    bool isViscosity { false };
+    bool isNumberOfDevices {false};
+    bool isDevices { false };
+    bool isOutputPath { false };
+    bool isPrefix { false };
+    bool isGridPath { false };
+    bool isPrintOutputFiles { false };
+    bool isGeometryValues { false };
+    bool isCalc2ndOrderMoments { false };
+    bool isCalc3rdOrderMoments { false };
+    bool isCalcHighOrderMoments { false };
+    bool isReadGeo { false };
+    bool isCalcMedian { false };
+    bool isCalcDragLift { false };
+    bool isCalcCp { false };
+    bool isWriteVeloASCII { false };
+    bool isCalcPlaneConc { false };
+    bool isConcFile { false };
+    bool isStreetVelocityFile { false };
+    bool isUseMeasurePoints { false };
+    bool isUseWale { false };
+    bool isUseInitNeq { false };
+    bool isSimulatePorousMedia { false };
+    bool isD3Qxx { false };
+    bool isTEnd { false };
+    bool isTOut { false };
+    bool isTStartOut { false };
+    bool isTimeCalcMedStart { false };
+    bool isTimeCalcMedEnd { false };
+    bool isPressInID { false };
+    bool isPressOutID { false };
+    bool isPressInZ { false };
+    bool isPressOutZ { false };
+    bool isDiffOn { false };
+    bool isDiffMod { false };
+    bool isDiffusivity { false };
+    bool isTemperatureInit { false };
+    bool isTemperatureBC { false };
+    // bool isViscosity { false };
+    bool isVelocity { false };
+    bool isViscosityRatio { false };
+    bool isVelocityRatio { false };
+    bool isDensityRatio { false };
+    bool isPressRatio { false };
+    bool isRealX { false };
+    bool isRealY { false };
+    bool isFactorPressBC { false };
+    bool isGeometryFileC { false };
+    bool isGeometryFileM { false };
+    bool isGeometryFileF { false };
+    bool isClockCycleForMP { false };
+    bool isTimestepForMP { false };
+    bool isForcingX { false };
+    bool isForcingY { false };
+    bool isForcingZ { false };
+    bool isQuadricLimiterP { false };
+    bool isQuadricLimiterM { false };
+    bool isQuadricLimiterD { false };
+    bool isCalcParticles { false };
+    bool isParticleBasicLevel { false };
+    bool isParticleInitLevel { false };
+    bool isNumberOfParticles { false };
+    bool isNeighborWSB { false };
+    bool isStartXHotWall { false };
+    bool isEndXHotWall { false };
+    bool isPossNeighborFilesX { false };
+    bool isPossNeighborFilesY { false };
+    bool isPossNeighborFilesZ { false };
+    bool isTimeDoCheckPoint { false };
+    bool isTimeDoRestart { false };
+    bool isDoCheckPoint { false };
+    bool isDoRestart { false };
+    bool isMaxLevel { false };
+    bool isGridX { false };
+    bool isGridY { false };
+    bool isGridZ { false };
+    bool isDistX { false };
+    bool isDistY { false };
+    bool isDistZ { false };
+    bool isNeedInterface { false };
+    bool isMainKernel { false };
+    bool isMultiKernelOn { false };
+    bool isMultiKernelLevel { false };
+    bool isMultiKernelName { false };
 };
 #endif
diff --git a/src/basics/basics/utilities/UbLogger.h b/src/basics/basics/utilities/UbLogger.h
index fc2b118715a0f4afc0251b8ed8e2373a7d488153..c0eb6bf243d52b6ed54a6b4f6cb7ac34ff24da64 100644
--- a/src/basics/basics/utilities/UbLogger.h
+++ b/src/basics/basics/utilities/UbLogger.h
@@ -217,8 +217,8 @@ inline std::string UbLogger<OutputPolicy>::logTimeString()
     }
 
     char result[100]   = { 0 };
-    static DWORD first = GetTickCount();
-    std::sprintf(result, "%s.%03ld", buffer, (long)(GetTickCount() - first) % 1000);
+    static DWORD first = GetTickCount64();
+    std::sprintf(result, "%s.%03ld", buffer, (long)(GetTickCount64() - first) % 1000);
     return result;
 }
 #else
diff --git a/src/basics/geometry3d/GbPolygon3D.cpp b/src/basics/geometry3d/GbPolygon3D.cpp
index 02138f75767420502d585842c7e05b9741df5b28..7905576ed5dae855350a333275cf5dcd1f6dbd3a 100644
--- a/src/basics/geometry3d/GbPolygon3D.cpp
+++ b/src/basics/geometry3d/GbPolygon3D.cpp
@@ -46,10 +46,13 @@ void GbPolygon3D::init()
 {
     x1s   = 0.0;
     x2s   = 0.0;
+    x3s   = 0.0;
     x1min = 0.0;
     x1max = 0.0;
     x2min = 0.0;
     x2max = 0.0;
+    x3min = 0.0;
+    x3max = 0.0;
     //		points   = NULL;
     consistent = false;
     ps         = NULL;
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index d60189b0747541e5e02fb38752df4581ee754b33..564cc2768469d0fdb83647126d2c94b6314d146b 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -204,6 +204,7 @@
 #include <CoProcessors/MPIIORestartCoProcessor.h>
 #include <CoProcessors/MicrophoneArrayCoProcessor.h>
 #include <WriteThixotropyQuantitiesCoProcessor.h>
+#include <WriteMultiphaseQuantitiesCoProcessor.h>
 
 #include <IntegrateValuesHelper.h>
 //#include <LBM/D3Q27CompactInterpolationProcessor.h>
@@ -239,6 +240,7 @@
 #include <LBM/RheologyK17LBMKernel.h>
 #include <LBM/RheologyPowellEyringModelLBMKernel.h>
 #include <LBM/MultiphaseCumulantLBMKernel.h>
+#include <LBM/MultiphaseScratchCumulantLBMKernel.h>
 
 
 
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..312ed01adf39ff6eb4aaf0965d8df6763ad3e8d1
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -0,0 +1,426 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file WriteMultiphaseQuantitiesCoProcessor.cpp
+//! \ingroup CoProcessors
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#include "WriteMultiphaseQuantitiesCoProcessor.h"
+#include "BCProcessor.h"
+#include "LBMKernel.h"
+#include <string>
+#include <vector>
+
+#include "BCArray3D.h"
+#include "Block3D.h"
+#include "Communicator.h"
+#include "DataSet3D.h"
+#include "Grid3D.h"
+#include "LBMUnitConverter.h"
+#include "UbScheduler.h"
+#include "basics/writer/WbWriterVtkXmlASCII.h"
+
+WriteMultiphaseQuantitiesCoProcessor::WriteMultiphaseQuantitiesCoProcessor() = default;
+//////////////////////////////////////////////////////////////////////////
+WriteMultiphaseQuantitiesCoProcessor::WriteMultiphaseQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
+                                                                             const std::string &path,
+                                                                             WbWriter *const writer,
+                                                                             SPtr<LBMUnitConverter> conv,
+                                                                             SPtr<Communicator> comm)
+        : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
+{
+    gridRank = comm->getProcessID();
+    minInitLevel = this->grid->getCoarsestInitializedLevel();
+    maxInitLevel = this->grid->getFinestInitializedLevel();
+
+    blockVector.resize(maxInitLevel + 1);
+
+    for (int level = minInitLevel; level <= maxInitLevel; level++)
+    {
+        grid->getBlocks(level, gridRank, true, blockVector[level]);
+    }
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::init()
+{}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::process(double step)
+{
+    if (scheduler->isDue(step))
+        collectData(step);
+
+    UBLOG(logDEBUG3, "WriteMultiphaseQuantitiesCoProcessor::update:" << step);
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::collectData(double step)
+{
+    int istep = static_cast<int>(step);
+
+    for (int level = minInitLevel; level <= maxInitLevel; level++)
+    {
+        for (SPtr<Block3D> block : blockVector[level])
+        {
+            if (block)
+            {
+                addDataMQ(block);
+            }
+        }
+    }
+
+    std::string pfilePath, partPath, subfolder, cfilePath;
+
+    subfolder = "mq" + UbSystem::toString(istep);
+    pfilePath = path + "/mq/" + subfolder;
+    cfilePath = path + "/mq/mq_collection";
+    partPath = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
+
+    std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
+    size_t found = partName.find_last_of("/");
+    std::string piece = partName.substr(found + 1);
+    piece = subfolder + "/" + piece;
+
+    std::vector<std::string> cellDataNames;
+    std::vector<std::string> pieces = comm->gather(piece);
+    if (comm->getProcessID() == comm->getRoot()) {
+        std::string pname =
+                WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
+        found = pname.find_last_of("/");
+        piece = pname.substr(found + 1);
+
+        std::vector<std::string> filenames;
+        filenames.push_back(piece);
+        if (step == CoProcessor::scheduler->getMinBegin())
+        {
+            WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
+        } else
+        {
+            WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
+        }
+        UBLOG(logINFO, "WriteMultiphaseQuantitiesCoProcessor step: " << istep);
+    }
+
+    clearData();
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::clearData()
+{
+    nodes.clear();
+    cells.clear();
+    datanames.clear();
+    data.clear();
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
+{
+    using namespace D3Q27System;
+    using namespace UbMath;
+
+    //double level   = (double)block->getLevel();
+
+    // Diese Daten werden geschrieben:
+    datanames.resize(0);
+    datanames.push_back("Phi");
+    datanames.push_back("Vx");
+    datanames.push_back("Vy");
+    datanames.push_back("Vz");
+    datanames.push_back("P1");
+
+    data.resize(datanames.size());
+
+    SPtr<LBMKernel> kernel                   = dynamicPointerCast<LBMKernel>(block->getKernel());
+    SPtr<BCArray3D> bcArray                  = kernel->getBCProcessor()->getBCArray();
+    SPtr<DistributionArray3D> distributionsF = kernel->getDataSet()->getFdistributions();
+    SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
+    SPtr<PhaseFieldArray3D> divU             = kernel->getDataSet()->getPhaseField();
+
+    LBMReal f[D3Q27System::ENDF + 1];
+    LBMReal phi[D3Q27System::ENDF + 1];
+    LBMReal vx1, vx2, vx3, rho, p1, beta, kappa;
+    LBMReal densityRatio = kernel->getDensityRatio();
+
+    kernel->getMultiphaseModelParameters(beta, kappa);
+    LBMReal phiL = kernel->getPhiL();
+    LBMReal phiH = kernel->getPhiH();
+
+    // knotennummerierung faengt immer bei 0 an!
+    int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+
+    if (block->getKernel()->getCompressible()) {
+        calcMacros = &D3Q27System::calcCompMacroscopicValues;
+    } else {
+        calcMacros = &D3Q27System::calcIncompMacroscopicValues;
+    }
+
+    // int minX1 = 0;
+    // int minX2 = 0;
+    // int minX3 = 0;
+
+    int maxX1 = (int)(distributionsF->getNX1());
+    int maxX2 = (int)(distributionsF->getNX2());
+    int maxX3 = (int)(distributionsF->getNX3());
+
+    int minX1 = 0;
+    int minX2 = 0;
+    int minX3 = 0;
+
+    // int maxX1 = (int)(distributions->getNX1());
+    // int maxX2 = (int)(distributions->getNX2());
+    // int maxX3 = (int)(distributions->getNX3());
+
+    // nummern vergeben und node vector erstellen + daten sammeln
+    CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
+    CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField1(
+        new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
+
+    for (int ix3 = minX3; ix3 < maxX3; ix3++) {
+        for (int ix2 = minX2; ix2 < maxX2; ix2++) {
+            for (int ix1 = minX1; ix1 < maxX1; ix1++) {
+                if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) {
+                    distributionsH->getDistribution(f, ix1, ix2, ix3);
+                    (*phaseField1)(ix1, ix2, ix3) =
+                        ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
+                        (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
+                         ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+                        ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+                }
+            }
+        }
+    }
+
+    maxX1 -= 2;
+    maxX2 -= 2;
+    maxX3 -= 2;
+
+    // maxX1 -= 1;
+    // maxX2 -= 1;
+    // maxX3 -= 1;
+
+    // D3Q27BoundaryConditionPtr bcPtr;
+    int nr = (int)nodes.size();
+    LBMReal dX1_phi;
+    LBMReal dX2_phi;
+    LBMReal dX3_phi;
+    LBMReal mu;
+
+    for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
+        for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
+            for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
+                if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) {
+                    int index                  = 0;
+                    nodeNumbers(ix1, ix2, ix3) = nr++;
+                    Vector3D worldCoordinates  = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+                    nodes.push_back(UbTupleFloat3(float(worldCoordinates[0]), float(worldCoordinates[1]),
+                                                  float(worldCoordinates[2])));
+
+                    phi[REST] = (*phaseField1)(ix1, ix2, ix3);
+
+                    if ((ix1 == 0) || (ix2 == 0) || (ix3 == 0)) {
+                        dX1_phi = 0.0;
+                        dX2_phi = 0.0;
+                        dX3_phi = 0.0;
+                        mu      = 0.0;
+                        // vx1 = 0.0;
+                        // vx2 = 0.0;
+                        // vx3 = 0.0;
+                    } else {
+                        phi[E]   = (*phaseField1)(ix1 + DX1[E], ix2 + DX2[E], ix3 + DX3[E]);
+                        phi[N]   = (*phaseField1)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
+                        phi[T]   = (*phaseField1)(ix1 + DX1[T], ix2 + DX2[T], ix3 + DX3[T]);
+                        phi[W]   = (*phaseField1)(ix1 + DX1[W], ix2 + DX2[W], ix3 + DX3[W]);
+                        phi[S]   = (*phaseField1)(ix1 + DX1[S], ix2 + DX2[S], ix3 + DX3[S]);
+                        phi[B]   = (*phaseField1)(ix1 + DX1[B], ix2 + DX2[B], ix3 + DX3[B]);
+                        phi[NE]  = (*phaseField1)(ix1 + DX1[NE], ix2 + DX2[NE], ix3 + DX3[NE]);
+                        phi[NW]  = (*phaseField1)(ix1 + DX1[NW], ix2 + DX2[NW], ix3 + DX3[NW]);
+                        phi[TE]  = (*phaseField1)(ix1 + DX1[TE], ix2 + DX2[TE], ix3 + DX3[TE]);
+                        phi[TW]  = (*phaseField1)(ix1 + DX1[TW], ix2 + DX2[TW], ix3 + DX3[TW]);
+                        phi[TN]  = (*phaseField1)(ix1 + DX1[TN], ix2 + DX2[TN], ix3 + DX3[TN]);
+                        phi[TS]  = (*phaseField1)(ix1 + DX1[TS], ix2 + DX2[TS], ix3 + DX3[TS]);
+                        phi[SW]  = (*phaseField1)(ix1 + DX1[SW], ix2 + DX2[SW], ix3 + DX3[SW]);
+                        phi[SE]  = (*phaseField1)(ix1 + DX1[SE], ix2 + DX2[SE], ix3 + DX3[SE]);
+                        phi[BW]  = (*phaseField1)(ix1 + DX1[BW], ix2 + DX2[BW], ix3 + DX3[BW]);
+                        phi[BE]  = (*phaseField1)(ix1 + DX1[BE], ix2 + DX2[BE], ix3 + DX3[BE]);
+                        phi[BS]  = (*phaseField1)(ix1 + DX1[BS], ix2 + DX2[BS], ix3 + DX3[BS]);
+                        phi[BN]  = (*phaseField1)(ix1 + DX1[BN], ix2 + DX2[BN], ix3 + DX3[BN]);
+                        phi[BSW] = (*phaseField1)(ix1 + DX1[BSW], ix2 + DX2[BSW], ix3 + DX3[BSW]);
+                        phi[BSE] = (*phaseField1)(ix1 + DX1[BSE], ix2 + DX2[BSE], ix3 + DX3[BSE]);
+                        phi[BNW] = (*phaseField1)(ix1 + DX1[BNW], ix2 + DX2[BNW], ix3 + DX3[BNW]);
+                        phi[BNE] = (*phaseField1)(ix1 + DX1[BNE], ix2 + DX2[BNE], ix3 + DX3[BNE]);
+                        phi[TNE] = (*phaseField1)(ix1 + DX1[TNE], ix2 + DX2[TNE], ix3 + DX3[TNE]);
+                        phi[TNW] = (*phaseField1)(ix1 + DX1[TNW], ix2 + DX2[TNW], ix3 + DX3[TNW]);
+                        phi[TSE] = (*phaseField1)(ix1 + DX1[TSE], ix2 + DX2[TSE], ix3 + DX3[TSE]);
+                        phi[TSW] = (*phaseField1)(ix1 + DX1[TSW], ix2 + DX2[TSW], ix3 + DX3[TSW]);
+                        dX1_phi  = 0.0 * gradX1_phi(phi);
+                        dX2_phi  = 0.0 * gradX2_phi(phi);
+                        dX3_phi  = 0.0 * gradX3_phi(phi);
+                        mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
+                    }
+
+                    distributionsF->getDistribution(f, ix1, ix2, ix3);
+                    //LBMReal dU = (*divU)(ix1, ix2, ix3);
+
+                    LBMReal rhoH = 1.0;
+                    LBMReal rhoL = 1.0 / densityRatio;
+                    // LBMReal rhoToPhi = (1.0 - 1.0/densityRatio);
+                    LBMReal rhoToPhi = (rhoH - rhoL) / (phiH - phiL);
+
+                    // rho = phi[ZERO] + (1.0 - phi[ZERO])*1.0/densityRatio;
+                    rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+
+                   vx1 =
+                        ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[BSE] - f[TNW]) + (f[BNE] - f[TSW]))) +
+                         (((f[BE] - f[TW]) + (f[TE] - f[BW])) + ((f[SE] - f[NW]) + (f[NE] - f[SW]))) + (f[E] - f[W])) /
+                            (rho * c1o3) +
+                        mu * dX1_phi / (2 * rho);
+
+                    vx2 =
+                        ((((f[TNE] - f[BSW]) + (f[BNW] - f[TSE])) + ((f[TNW] - f[BSE]) + (f[BNE] - f[TSW]))) +
+                         (((f[BN] - f[TS]) + (f[TN] - f[BS])) + ((f[NW] - f[SE]) + (f[NE] - f[SW]))) + (f[N] - f[S])) /
+                            (rho * c1o3) +
+                        mu * dX2_phi / (2 * rho);
+
+                    vx3 =
+                        ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[TNW] - f[BSE]) + (f[TSW] - f[BNE]))) +
+                         (((f[TS] - f[BN]) + (f[TN] - f[BS])) + ((f[TW] - f[BE]) + (f[TE] - f[BW]))) + (f[T] - f[B])) /
+                            (rho * c1o3) +
+                        mu * dX3_phi / (2 * rho);
+
+                    p1 = (((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
+                          (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
+                           ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+                          ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST]) +
+                         (vx1 * rhoToPhi * dX1_phi * c1o3 + vx2 * rhoToPhi * dX2_phi * c1o3 +
+                          vx3 * rhoToPhi * dX3_phi * c1o3) /
+                             2.0;
+                    p1 = rho * p1 * c1o3;
+
+                    // calcMacros(f,rho,vx1,vx2,vx3);
+                    // tempfield(ix1,ix2,ix3)= ; Added by HS
+                    // double press = D3Q27System::calcPress(f,rho,vx1,vx2,vx3);
+
+					if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
+                        UB_THROW(UbException(
+                            UB_EXARGS, "vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                           block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+                                           UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+                    // vx1=999.0;
+                    if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
+                        UB_THROW(UbException(
+                            UB_EXARGS, "vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                           block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+                                           UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+                    // vx2=999.0;
+                    if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
+                        UB_THROW(UbException(
+                            UB_EXARGS, "vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                           block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+                                           UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+
+                    if (UbMath::isNaN(phi[REST]) || UbMath::isInfinity(phi[REST]))
+                        UB_THROW(UbException(
+                            UB_EXARGS, "phi is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                           block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+                                           UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+
+                    if (UbMath::isNaN(p1) || UbMath::isInfinity(p1))
+                        UB_THROW( UbException(UB_EXARGS,"p1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
+                        ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+
+                    data[index++].push_back(phi[REST]);
+                    data[index++].push_back(vx1);
+                    data[index++].push_back(vx2);
+                    data[index++].push_back(vx3);
+                    data[index++].push_back(p1);
+                }
+            }
+        }
+    }
+    maxX1 -= 1;
+    maxX2 -= 1;
+    maxX3 -= 1;
+    // cell vector erstellen
+    for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
+        for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
+            for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
+                if ((SWB = nodeNumbers(ix1, ix2, ix3)) >= 0 && (SEB = nodeNumbers(ix1 + 1, ix2, ix3)) >= 0 &&
+                    (NEB = nodeNumbers(ix1 + 1, ix2 + 1, ix3)) >= 0 && (NWB = nodeNumbers(ix1, ix2 + 1, ix3)) >= 0 &&
+                    (SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
+                    (NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
+                    (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
+                    cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB,
+                                                (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET,
+                                                (unsigned int)NET, (unsigned int)NWT));
+                }
+            }
+        }
+    }
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX1_phi(const LBMReal *const &h)
+{
+    using namespace D3Q27System;
+    LBMReal sum = 0.0;
+    for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+        sum += WEIGTH[k] * DX1[k] * h[k];
+    }
+    return 3.0 * sum;
+}
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX2_phi(const LBMReal *const &h)
+{
+    using namespace D3Q27System;
+    LBMReal sum = 0.0;
+    for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+        sum += WEIGTH[k] * DX2[k] * h[k];
+    }
+    return 3.0 * sum;
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX3_phi(const LBMReal *const &h)
+{
+    using namespace D3Q27System;
+    LBMReal sum = 0.0;
+    for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+        sum += WEIGTH[k] * DX3[k] * h[k];
+    }
+    return 3.0 * sum;
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::nabla2_phi(const LBMReal *const &h)
+{
+    using namespace D3Q27System;
+    LBMReal sum = 0.0;
+    for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+        sum += WEIGTH[k] * (h[k] - h[REST]);
+    }
+    return 6.0 * sum;
+}
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4504c1dfccbad377e0bad4bc1aab51989abaff4
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h
@@ -0,0 +1,104 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file WriteMultiphaseQuantitiesCoProcessor.h
+//! \ingroup CoProcessors
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#ifndef WriteMultiphaseQuantitiesCoProcessor_H
+#define WriteMultiphaseQuantitiesCoProcessor_H
+
+#include <PointerDefinitions.h>
+#include <string>
+#include <vector>
+
+#include "CoProcessor.h"
+#include "LBMSystem.h"
+#include "UbTuple.h"
+
+class Communicator;
+class Grid3D;
+class UbScheduler;
+class LBMUnitConverter;
+class WbWriter;
+class Block3D;
+
+//! \brief A class writes macroscopic quantities information to a VTK-file
+class WriteMultiphaseQuantitiesCoProcessor : public CoProcessor
+{
+public:
+    WriteMultiphaseQuantitiesCoProcessor();
+    //! \brief Construct WriteMultiphaseQuantitiesCoProcessor object
+    //! \pre The Grid3D and UbScheduler objects must exist
+    //! \param grid is observable Grid3D object
+    //! \param s is UbScheduler object for scheduling of observer
+    //! \param path is path of folder for output
+    //! \param writer is WbWriter object
+    //! \param conv is LBMUnitConverter object
+    //! \param comm is Communicator object
+    WriteMultiphaseQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path,
+                                          WbWriter *const writer, SPtr<LBMUnitConverter> conv, SPtr<Communicator> comm);
+    ~WriteMultiphaseQuantitiesCoProcessor() override = default;
+
+    void process(double step) override;
+
+protected:
+    //! Collect data for VTK-file
+    //! \param step is a time step
+    void collectData(double step);
+    //! Collect data for VTK-file
+    //! \param block is a time step
+    void addDataMQ(SPtr<Block3D> block);
+    void clearData();
+
+private:
+    void init();
+    std::vector<UbTupleFloat3> nodes;
+    std::vector<UbTupleUInt8> cells;
+    std::vector<std::string> datanames;
+    std::vector<std::vector<double>> data;
+    std::string path;
+    WbWriter *writer;
+    SPtr<LBMUnitConverter> conv;
+    std::vector<std::vector<SPtr<Block3D>>> blockVector;
+    int minInitLevel;
+    int maxInitLevel;
+    int gridRank;
+    SPtr<Communicator> comm;
+
+    LBMReal gradX1_phi(const LBMReal *const &);
+    LBMReal gradX2_phi(const LBMReal *const &);
+    LBMReal gradX3_phi(const LBMReal *const &);
+    LBMReal nabla2_phi(const LBMReal *const &);
+
+    using CalcMacrosFct = void (*)(const LBMReal *const &, LBMReal &, LBMReal &, LBMReal &, LBMReal &);
+    CalcMacrosFct calcMacros;
+};
+
+#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index 1bfa01272911c6da84208a9f8484ea6f2be5b78c..f8119ea3557cc528407c560f415c0d7115927a60 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -207,7 +207,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
         }
 
         LBMReal collFactorM;
-        LBMReal forcingTerm[D3Q27System::ENDF + 1];
+        //LBMReal forcingTerm[D3Q27System::ENDF + 1];
 
         for (int x3 = minX3; x3 < maxX3; x3++) {
             for (int x2 = minX2; x2 < maxX2; x2++) {
@@ -286,7 +286,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
 
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        //LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
                         LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);