From 6a96317d60ed0645b347fcbc2c49fd7ab0006e6b Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Fri, 17 Sep 2021 13:04:25 +0200
Subject: [PATCH] fix bug and remove redundant initialisation in
 MultiphaseSetKernelBlockVisitor / remove unused parameter in
 MultiphaseInitDistributionsBlockVisitor

---
 apps/cpu/Multiphase/Multiphase.cpp            |  2 +-
 .../cpu/MultiphaseDropletTest/DropletTest.cfg | 33 ++++++------
 apps/cpu/MultiphaseDropletTest/droplet.cpp    | 45 ++++++++++++----
 src/cpu/VirtualFluids.h                       |  2 +-
 ...eTwoPhaseFieldsPressureFilterLBMKernel.cpp | 19 +++----
 ...woPhaseFieldsVelocityCumulantLBMKernel.cpp |  2 +-
 .../Parallel/MPICommunicator.h                | 21 +++++---
 ...ultiphaseInitDistributionsBlockVisitor.cpp |  7 +--
 .../MultiphaseInitDistributionsBlockVisitor.h |  4 +-
 .../MultiphaseSetKernelBlockVisitor.cpp       |  4 +-
 ...cityFormInitDistributionsBlockVisitor.cpp} | 53 ++++++++-----------
 ...locityFormInitDistributionsBlockVisitor.h} | 16 +++---
 12 files changed, 109 insertions(+), 99 deletions(-)
 rename src/cpu/VirtualFluidsCore/Visitors/{MultiphaseInitDistributionsBlockVisitorVelocity.cpp => MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp} (84%)
 rename src/cpu/VirtualFluidsCore/Visitors/{MultiphaseInitDistributionsBlockVisitorVelocity.h => MultiphaseVelocityFormInitDistributionsBlockVisitor.h} (86%)

diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index 187e655eb..4c5904a21 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -363,7 +363,7 @@ void run(string configname)
             mu::Parser fct1;
             fct1.SetExpr("phiL");
             fct1.DefineConst("phiL", phiL);
-            MultiphaseInitDistributionsBlockVisitor initVisitor(densityRatio, interfaceThickness, radius);
+            MultiphaseInitDistributionsBlockVisitor initVisitor(interfaceThickness);
             initVisitor.setPhi(fct1);
             grid->accept(initVisitor);
 
diff --git a/apps/cpu/MultiphaseDropletTest/DropletTest.cfg b/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
index 243e04299..78b53897a 100644
--- a/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
+++ b/apps/cpu/MultiphaseDropletTest/DropletTest.cfg
@@ -1,36 +1,33 @@
 #pathname = d:/temp/MultiphaseDropletTest
-pathname = E:/Multiphase/DropletTestSigma2
+pathname = E:/Multiphase/DropletTest
 
-numOfThreads = 1
+numOfThreads = 4
 availMem = 10e9
 
 #Grid
 
-boundingBox = 0 128 0 64 0 64
-#boundingBox = 0 8 0 8 0 8
-blocknx = 8 8 8
-#blocknx = 64 32 32
-#blocknx= 128 64 64
+boundingBox = 0 50 0 150 0 50
+blocknx = 10 10 10
 
 dx = 1
 refineLevel = 0
 
 #Simulation
-uLB = 0.01#0.005#0.005 
+uLB = 0#0.01#0.005#0.005 
 Re = 10
-nuL =1e-3# 1.0e-5 #!1e-2
-nuG =1e-6# 1.16e-4 #!1e-2
-densityRatio = 1000
-sigma =0.0001# 1e-5 #4.66e-3 #surface tension 1e-4 ./. 1e-5
+nuL = 1e-3 #1e-5# 1.0e-5 #!1e-2
+nuG = 1e-6 # 1e-8 # 1.16e-4 #!1e-2
+densityRatio = 10
+sigma = 1.111111111111111e-07 #1e-10 #1e-6  # 1e-5 #4.66e-3 #surface tension 1e-4 ./. 1e-5
 interfaceThickness = 5
-radius = 8#16
+radius = 5 #16
 contactAngle = 110.0
-gravity = 0.0
-#gravity = 5.04e-6
+#gravity = 0.0
+gravity = -1.111111111111111e-07 #-5.04e-6
 phi_L = 0.0
 phi_H = 1.0
 Phase-field Relaxation = 0.6
-Mobility = 0.02 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
+Mobility = 0.1 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
 
 
 logToFile = false
@@ -41,5 +38,5 @@ restartStep = 100000
 cpStart = 100000
 cpStep = 100000
 
-outTime = 1
-endTime = 10000
\ No newline at end of file
+outTime = 100
+endTime = 1000000
\ No newline at end of file
diff --git a/apps/cpu/MultiphaseDropletTest/droplet.cpp b/apps/cpu/MultiphaseDropletTest/droplet.cpp
index a6892cad1..e022ce0c3 100644
--- a/apps/cpu/MultiphaseDropletTest/droplet.cpp
+++ b/apps/cpu/MultiphaseDropletTest/droplet.cpp
@@ -72,6 +72,8 @@ void run(string configname)
         LBMReal rhoLB = 0.0;
         LBMReal nuLB  = nuL; //(uLB*dLB) / Re;
 
+        LBMReal D  = 2.0*radius;
+
         SPtr<LBMUnitConverter> conv(new LBMUnitConverter());
 
         //const int baseLevel = 0;
@@ -82,9 +84,10 @@ void run(string configname)
         //kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
         kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
 
+
         kernel->setWithForcing(true);
-        kernel->setForcingX1(gr);
-        kernel->setForcingX2(0.0);
+        kernel->setForcingX1(0.0);
+        kernel->setForcingX2(gr);
         kernel->setForcingX3(0.0);
 
         kernel->setPhiL(phiL);
@@ -92,16 +95,28 @@ void run(string configname)
         kernel->setPhaseFieldRelaxation(tauH);
         kernel->setMobility(mob);
 
+        kernel->setCollisionFactorMultiphase(nuL, nuG);
+        kernel->setDensityRatio(densityRatio);
+        kernel->setMultiphaseModelParameters(beta, kappa);
+        kernel->setContactAngle(theta);
+
         SPtr<BCProcessor> bcProc(new BCProcessor());
         // BCProcessorPtr bcProc(new ThinWallBCProcessor());
 
         kernel->setBCProcessor(bcProc);
 
+        SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+        noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseNoSlipBCAlgorithm()));
+        //////////////////////////////////////////////////////////////////////////////////
+        // BC visitor
+        MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
+        bcVisitor.addBC(noSlipBCAdapter);
+
         SPtr<Grid3D> grid(new Grid3D(comm));
         grid->setDeltaX(dx);
         grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
         grid->setPeriodicX1(true);
-        grid->setPeriodicX2(true);
+        grid->setPeriodicX2(false);
         grid->setPeriodicX3(true);
         grid->setGhostLayerWidth(2);
 
@@ -148,6 +163,14 @@ void run(string configname)
             GenBlocksGridVisitor genBlocks(gridCube);
             grid->accept(genBlocks);
 
+            double dx2 = 2.0 * dx;
+            GbCuboid3DPtr wallYmin(new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_minX2, g_maxX3 + dx2));
+            GbSystem3D::writeGeoObject(wallYmin.get(), pathname + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
+            GbCuboid3DPtr wallYmax(new GbCuboid3D(g_minX1 - dx2, g_maxX2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
+            GbSystem3D::writeGeoObject(wallYmax.get(), pathname + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
+
+            SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+            SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
  
             SPtr<WriteBlocksCoProcessor> ppblocks(new WriteBlocksCoProcessor(
                 grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
@@ -155,6 +178,8 @@ void run(string configname)
             SPtr<Grid3DVisitor> metisVisitor(
                 new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
             InteractorsHelper intHelper(grid, metisVisitor, true);
+            intHelper.addInteractor(wallYminInt);
+            intHelper.addInteractor(wallYmaxInt);
             intHelper.selectBlocks();
 
             ppblocks->process(0);
@@ -199,12 +224,10 @@ void run(string configname)
 
             intHelper.setBC();
 
-            //grid->accept(bcVisitor);
-
             // initialization of distributions
-            LBMReal x1c = (g_maxX1 - g_minX1-1)/2;
-            LBMReal x2c = (g_maxX2 - g_minX2-1)/2;
-            LBMReal x3c = (g_maxX3 - g_minX3-1)/2;
+            LBMReal x1c = 2.5 * D; // (g_maxX1 - g_minX1-1)/2; //
+            LBMReal x2c = 12.5 * D; //(g_maxX2 - g_minX2-1)/2;
+            LBMReal x3c = 2.5 * D; //(g_maxX3 - g_minX3-1)/2;
             mu::Parser fct1;
             fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
             fct1.DefineConst("x1c", x1c);
@@ -222,7 +245,7 @@ void run(string configname)
             fct2.DefineConst("radius", radius);
             fct2.DefineConst("interfaceThickness", interfaceThickness);
 
-            MultiphaseInitDistributionsBlockVisitorVelocity initVisitor(densityRatio, interfaceThickness, radius);
+            MultiphaseInitDistributionsBlockVisitor initVisitor(densityRatio);
             initVisitor.setPhi(fct1);
             initVisitor.setVx1(fct2);
             grid->accept(initVisitor);
@@ -232,7 +255,7 @@ void run(string configname)
                 SPtr<UbScheduler> geoSch(new UbScheduler(1));
                 SPtr<WriteBoundaryConditionsCoProcessor> ppgeo(new WriteBoundaryConditionsCoProcessor(
                     grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
-                //ppgeo->process(0);
+                ppgeo->process(0);
                 ppgeo.reset();
             }
 
@@ -258,6 +281,8 @@ void run(string configname)
                 UBLOG(logINFO, "Restart - end");
         }
 
+        grid->accept(bcVisitor);
+
         //TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
         //grid->accept(setConnsVisitor);
 
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index 315de6dd7..7edd56779 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -341,7 +341,7 @@
 #include <Visitors/MultiphaseSetKernelBlockVisitor.h>
 #include <Visitors/MultiphaseBoundaryConditionsBlockVisitor.h>
 #include <Visitors/MultiphaseInitDistributionsBlockVisitor.h>
-#include <Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.h>
+#include <Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.h>
 #include <Visitors/SetInterpolationConnectorsBlockVisitor.h>
 
 #include <RefineAroundGbObjectHelper.h>
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
index f64f2719c..6ffbdb498 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
@@ -139,7 +139,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
     forcingX2 = 0.0;
     forcingX3 = 0.0;
 
-	LBMReal oneOverInterfaceScale = 1.5;
+	LBMReal oneOverInterfaceScale = 1.0;//1.5;
     /////////////////////////////////////
 
     localDistributionsF    = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
@@ -215,7 +215,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
                         LBMReal mfcca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
                         LBMReal mfbbb = (*this->zeroDistributionsH1)(x1, x2, x3);
-                        (*phaseField)(x1, x2, x3) = (((mfaaa + mfccc) + (mfaca + mfcac)) + ((mfaac + mfcca)  + (mfcaa + mfacc))  ) +
+						LBMReal testPhi1 = (*phaseField)(x1, x2, x3) = (((mfaaa + mfccc) + (mfaca + mfcac)) + ((mfaac + mfcca)  + (mfcaa + mfacc))  ) +
                                                     (((mfaab + mfacb) + (mfcab + mfccb)) + ((mfaba + mfabc) + (mfcba + mfcbc)) +
                                                     ((mfbaa + mfbac) + (mfbca + mfbcc))) + ((mfabb + mfcbb) +
                                                     (mfbab + mfbcb) + (mfbba + mfbbc)) + mfbbb;
@@ -248,7 +248,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
                         mfcca = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BNE, x1, x2, x3p);
 
                         mfbbb = (*this->zeroDistributionsH2)(x1, x2, x3);
-                        (*phaseField2)(x1, x2, x3) =
+						LBMReal testPhi2 = (*phaseField2)(x1, x2, x3) =
                             (((mfaaa + mfccc) + (mfaca + mfcac)) + ((mfaac + mfcca) + (mfcaa + mfacc))) +
                             (((mfaab + mfacb) + (mfcab + mfccb)) + ((mfaba + mfabc) + (mfcba + mfcbc)) +
                              ((mfbaa + mfbac) + (mfbca + mfbcc))) +
@@ -305,8 +305,9 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 
 						(*pressure)(x1, x2, x3) = (*pressure)(x1, x2, x3) + rho * c1o3 * drho;
 
+						//(*pressure)(x1, x2, x3) = (((*phaseField)(x1, x2, x3)) + ((*phaseField2)(x1, x2, x3)) - c1) * c1o3;
 						////!!!!!! relplace by pointer swap!
-						(*pressureOld)(x1, x2, x3) = (*pressure)(x1, x2, x3);
+						//(*pressureOld)(x1, x2, x3) = (*pressure)(x1, x2, x3);
                     }
                 }
             }
@@ -603,9 +604,9 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 			   LBMReal rhoRef=c1;
 
 			  //LBMReal 
-			//    LBMReal drho = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
-			// 	   + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
-			// 	   + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
+			   LBMReal drho = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
+				   + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
+				   + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
 
 			   LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
 				   (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
@@ -701,7 +702,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 
 			   //Viscosity increase by pressure gradient
 			   LBMReal errPhi = (((1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale)- denom);
-			   LBMReal limVis = 0.0000001;//0.01;
+			   LBMReal limVis = 0.0000001*10;//0.01;
 			  // collFactorM =collFactorM/(c1+limVis*(errPhi*errPhi)*collFactorM);
 			  // collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
 			   errPhi = errPhi * errPhi* errPhi * errPhi * errPhi * errPhi;
@@ -1250,7 +1251,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 
 			   LBMReal OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
 			   LBMReal OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
-			//    LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
+			   LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
 			   LBMReal A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
 			   //FIXME:  warning C4459: declaration of 'B' hides global declaration (message : see declaration of 'D3Q27System::B' )
 			   LBMReal BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
index faa9bd438..294568b4f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
@@ -714,7 +714,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::calculate(int step)
 
 
 			   ///surface tension force
-			   vvx += mu * dX1_phi*c1o2;
+			   vvx += mu * dX1_phi * c1o2;
 			   vvy += mu * dX2_phi * c1o2 ;
 			   vvz += mu * dX3_phi * c1o2;
 
diff --git a/src/cpu/VirtualFluidsCore/Parallel/MPICommunicator.h b/src/cpu/VirtualFluidsCore/Parallel/MPICommunicator.h
index 3c25cc78f..19ad23a53 100644
--- a/src/cpu/VirtualFluidsCore/Parallel/MPICommunicator.h
+++ b/src/cpu/VirtualFluidsCore/Parallel/MPICommunicator.h
@@ -139,16 +139,23 @@ void MPICommunicator::allGather(std::vector<T> &svalues, std::vector<T> &rvalues
 
     rvalues.resize(displs[numprocs - 1] + rcounts[numprocs - 1]);
 
-    if (rvalues.size() == 0) {
-        rvalues.resize(1);
-        rvalues[0] = 999;
+    T* sval = NULL;
+    T* rval = NULL;
+
+    if (svalues.size() > 0) {
+        //svalues.resize(1);
+        //svalues[0] = 999;
+        sval = &svalues[0];
     }
-    if (scount == 0) {
-        svalues.resize(1);
-        svalues[0] = 999;
+
+    if (rvalues.size() > 0) {
+        //rvalues.resize(1);
+        //rvalues[0] = 999;
+        rval = &rvalues[0];
     }
 
-    MPI_Allgatherv(&svalues[0], scount, mpiDataType, &rvalues[0], &rcounts[0], &displs[0], mpiDataType, comm);
+    //MPI_Allgatherv(&svalues[0], scount, mpiDataType, &rvalues[0], &rcounts[0], &displs[0], mpiDataType, comm);
+    MPI_Allgatherv(sval, scount, mpiDataType, rval, &rcounts[0], &displs[0], mpiDataType, comm);
 }
 //////////////////////////////////////////////////////////////////////////
 template <class T>
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
index 124e342d6..a59e8d476 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
@@ -50,12 +50,9 @@ MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor
 	this->setRho(0.0);
 }
 //////////////////////////////////////////////////////////////////////////
-MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor( LBMReal densityRatio, LBMReal intThickness, LBMReal radius, LBMReal vx1, LBMReal vx2, LBMReal vx3)
-	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), densityRatio(densityRatio) /*, intThickness(intThickness), radius(radius) */
+MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor( LBMReal densityRatio, LBMReal vx1, LBMReal vx2, LBMReal vx3)
+	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), densityRatio(densityRatio) 
 {
-    (void) intThickness;
-    (void) radius;
-
 	this->setVx1(vx1);
 	this->setVx2(vx2);
 	this->setVx3(vx3);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.h
index f94167040..e351b901c 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.h
@@ -56,7 +56,7 @@ public:
 	//! \param vx1 - velocity in x
 	//! \param vx2 - velocity in y
 	//! \param vx3 - velocity in z
-	MultiphaseInitDistributionsBlockVisitor( LBMReal densityRatio, LBMReal intThickness, LBMReal radius, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0);
+	MultiphaseInitDistributionsBlockVisitor( LBMReal densityRatio, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0);
 	//////////////////////////////////////////////////////////////////////////
 	//automatic vars are: x1,x2, x3
 	//ussage example: setVx1("x1*0.01+x2*0.003")
@@ -95,8 +95,6 @@ private:
 
 	LBMReal nu;
 	LBMReal densityRatio;
-	//LBMReal intThickness;
-	//LBMReal radius;
 };
 
 #endif //D3Q27INITDISTRIBUTIONSPATCHVISITOR_H
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
index e42cdb785..11863b691 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
@@ -35,8 +35,8 @@ void MultiphaseSetKernelBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> blo
 	{
 		LBMReal collFactorL = LBMSystem::calcCollisionFactor(nuL, block->getLevel());
 		LBMReal collFactorG = LBMSystem::calcCollisionFactor(nuG, block->getLevel());
-
-		//kernel->setCollisionFactorMultiphase(collFactorL, collFactorG);
+		kernel->setCollisionFactorMultiphase(collFactorL, collFactorG);
+		
 		//kernel->setDensityRatio(densityRatio);
 		//kernel->setMultiphaseModelParameters(beta, kappa);
 		//kernel->setContactAngle(contactAngle);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
similarity index 84%
rename from src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.cpp
rename to src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
index bd08dada8..8a201fa3e 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
@@ -26,12 +26,12 @@
 //  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 MultiphaseInitDistributionsBlockVisitorVelocity.cpp
+//! \file MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
 //! \ingroup Visitors
 //! \author Hesameddin Safari
 //=======================================================================================
 
-#include "MultiphaseInitDistributionsBlockVisitorVelocity.h"
+#include "MultiphaseVelocityFormInitDistributionsBlockVisitor.h"
 #include "BCArray3D.h"
 #include "BCProcessor.h"
 #include "Block3D.h"
@@ -41,7 +41,7 @@
 #include "Grid3DSystem.h"
 #include "LBMKernel.h"
 
-MultiphaseInitDistributionsBlockVisitorVelocity::MultiphaseInitDistributionsBlockVisitorVelocity() 
+MultiphaseVelocityFormInitDistributionsBlockVisitor::MultiphaseVelocityFormInitDistributionsBlockVisitor() 
 	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
 {
 	this->setVx1(0.0);
@@ -50,108 +50,97 @@ MultiphaseInitDistributionsBlockVisitorVelocity::MultiphaseInitDistributionsBloc
 	this->setRho(0.0);
 }
 //////////////////////////////////////////////////////////////////////////
-MultiphaseInitDistributionsBlockVisitorVelocity::MultiphaseInitDistributionsBlockVisitorVelocity( /*LBMReal densityRatio,*/ LBMReal intThickness, LBMReal radius, LBMReal vx1, LBMReal vx2, LBMReal vx3)
-	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL)/*, densityRatio(densityRatio), intThickness(intThickness), radius(radius) */
-{
-    (void) intThickness;
-    (void) radius;
-
-	this->setVx1(vx1);
-	this->setVx2(vx2);
-	this->setVx3(vx3);
-}
-//////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx1( const mu::Parser& parser)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx1( const mu::Parser& parser)  
 { 
 	this->checkFunction(parser); 
 	this->muVx1 = parser;  
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx2( const mu::Parser& parser)
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx2( const mu::Parser& parser)
 { 
 	this->checkFunction(parser); 
 	this->muVx2 = parser;  
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx3( const mu::Parser& parser)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx3( const mu::Parser& parser)  
 { 
 	this->checkFunction(parser); 
 	this->muVx3 = parser;  
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setRho( const mu::Parser& parser)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setRho( const mu::Parser& parser)  
 { 
 	this->checkFunction(parser); 
 	this->muRho = parser;  
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setPhi( const mu::Parser& parser)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setPhi( const mu::Parser& parser)  
 { 
 	this->checkFunction(parser); 
 	this->muPhi = parser;  
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx1( const std::string& muParserString)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx1( const std::string& muParserString)  
 { 
 	this->muVx1.SetExpr(muParserString); 
 	this->checkFunction(muVx1); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx2( const std::string& muParserString) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx2( const std::string& muParserString) 
 { 
 	this->muVx2.SetExpr(muParserString); 
 	this->checkFunction(muVx2); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx3( const std::string& muParserString)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx3( const std::string& muParserString)  
 { 
 	this->muVx3.SetExpr(muParserString); 
 	this->checkFunction(muVx3); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setRho( const std::string& muParserString)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setRho( const std::string& muParserString)  
 { 
 	this->muRho.SetExpr(muParserString); 
 	this->checkFunction(muRho); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setPhi( const std::string& muParserString)  
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setPhi( const std::string& muParserString)  
 { 
 	this->muPhi.SetExpr(muParserString); 
 	this->checkFunction(muPhi); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx1( LBMReal vx1 ) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx1( LBMReal vx1 ) 
 { 
 	this->muVx1.SetExpr( UbSystem::toString(vx1,D3Q27RealLim::digits10) );  
 	this->checkFunction(muVx1); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx2( LBMReal vx2 ) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx2( LBMReal vx2 ) 
 { 
 	this->muVx2.SetExpr( UbSystem::toString(vx2,D3Q27RealLim::digits10) );  
 	this->checkFunction(muVx2); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setVx3( LBMReal vx3 ) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setVx3( LBMReal vx3 ) 
 { 
 	this->muVx3.SetExpr( UbSystem::toString(vx3,D3Q27RealLim::digits10) );  
 	this->checkFunction(muVx3); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setRho( LBMReal rho ) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setRho( LBMReal rho ) 
 { 
 	this->muRho.SetExpr( UbSystem::toString(rho,D3Q27RealLim::digits10) );  
 	this->checkFunction(muRho); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setPhi( LBMReal phi ) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setPhi( LBMReal phi ) 
 { 
 	this->muPhi.SetExpr( UbSystem::toString(phi,D3Q27RealLim::digits10) );  
 	this->checkFunction(muPhi); 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::visit(const SPtr<Grid3D> grid, SPtr<Block3D> block) 
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> block) 
 {
 	using namespace D3Q27System;
 
@@ -336,7 +325,7 @@ void MultiphaseInitDistributionsBlockVisitorVelocity::visit(const SPtr<Grid3D> g
 
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::checkFunction(mu::Parser fct)
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::checkFunction(mu::Parser fct)
 {
 	double x1=1.0,x2=1.0,x3=1.0;
 	fct.DefineVar("x1",&x1); 
@@ -355,7 +344,7 @@ void MultiphaseInitDistributionsBlockVisitorVelocity::checkFunction(mu::Parser f
 	}
 }
 //////////////////////////////////////////////////////////////////////////
-void MultiphaseInitDistributionsBlockVisitorVelocity::setNu( LBMReal nu )
+void MultiphaseVelocityFormInitDistributionsBlockVisitor::setNu( LBMReal nu )
 {
 	this->nu = nu;
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.h b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.h
similarity index 86%
rename from src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.h
rename to src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.h
index 5426817d5..72841a7c5 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitorVelocity.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.h
@@ -26,13 +26,13 @@
 //  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 MultiphaseInitDistributionsBlockVisitorVelocity.h
+//! \file MultiphaseVelocityFormInitDistributionsBlockVisitor.h
 //! \ingroup Visitors
-//! \author Hesameddin Safari
+//! \author Hesameddin Safari, Martin Geier, Konstantin Kutscher
 //=======================================================================================
 
-#ifndef MultiphaseInitDistributionsBlockVisitorVelocity_H
-#define MultiphaseInitDistributionsBlockVisitorVelocity_H
+#ifndef MultiphaseVelocityFormInitDistributionsBlockVisitor_H
+#define MultiphaseVelocityFormInitDistributionsBlockVisitor_H
 
 #include "Block3DVisitor.h"
 #include "D3Q27System.h"
@@ -42,13 +42,13 @@
 
 
 
-class MultiphaseInitDistributionsBlockVisitorVelocity : public Block3DVisitor
+class MultiphaseVelocityFormInitDistributionsBlockVisitor : public Block3DVisitor
 {
 public:
 	typedef std::numeric_limits<LBMReal> D3Q27RealLim;
 
 public:
-	MultiphaseInitDistributionsBlockVisitorVelocity();
+	MultiphaseVelocityFormInitDistributionsBlockVisitor();
 	//D3Q27ETInitDistributionsBlockVisitor(LBMReal rho, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0);
 	//! Constructor
 	//! \param nu - viscosity
@@ -56,7 +56,6 @@ public:
 	//! \param vx1 - velocity in x
 	//! \param vx2 - velocity in y
 	//! \param vx3 - velocity in z
-	MultiphaseInitDistributionsBlockVisitorVelocity( /*LBMReal densityRatio, */LBMReal intThickness, LBMReal radius, LBMReal vx1=0.0, LBMReal vx2=0.0, LBMReal vx3=0.0);
 	//////////////////////////////////////////////////////////////////////////
 	//automatic vars are: x1,x2, x3
 	//ussage example: setVx1("x1*0.01+x2*0.003")
@@ -94,9 +93,6 @@ private:
 	mu::Parser muPhi;
 
 	LBMReal nu;
-	//LBMReal densityRatio;
-	//LBMReal intThickness;
-	//LBMReal radius;
 };
 
 #endif //D3Q27INITDISTRIBUTIONSPATCHVISITOR_H
-- 
GitLab