From 9e1c774aebe6d7af165be2bd52d690fa8d506bb7 Mon Sep 17 00:00:00 2001
From: alena <akaranchuk@list.ru>
Date: Mon, 28 Mar 2022 12:01:17 +0200
Subject: [PATCH] fix some problems: deadlock und nan reasons

---
 apps/cpu/RisingBubble2D/RisingBubble2D.cpp    | 101 +++++++-----------
 .../LBM/MultiphasePressureFilterLBMKernel.cpp |  75 ++++++++++++-
 2 files changed, 109 insertions(+), 67 deletions(-)

diff --git a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
index 0182f569c..7d8361e00 100644
--- a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
+++ b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
@@ -46,7 +46,7 @@ void run(string configname)
         double cpStart     = config.getValue<double>("cpStart");
         double cpStep      = config.getValue<double>("cpStep");
         bool newStart      = config.getValue<bool>("newStart");
-        //double rStep = config.getValue<double>("rStep");
+        double rStep = config.getValue<double>("rStep");
 
         std::shared_ptr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
         int myid                = comm->getProcessID();
@@ -85,7 +85,7 @@ void run(string configname)
 //         }    
 //#endif
 
-        //Sleep(30000);
+        //Sleep(20000);
 
         // LBMReal dLB = 0; // = length[1] / dx;
         LBMReal rhoLB = 0.0;
@@ -95,18 +95,18 @@ void run(string configname)
         LBMReal D  = 2.0*radius;
 
         //density retio
-        //LBMReal r_rho = densityRatio;
+        LBMReal r_rho = densityRatio;
 
         //density of heavy fluid
         LBMReal rho_h = 1.0;
         //density of light fluid
-        //LBMReal rho_l = rho_h / r_rho;
+        LBMReal rho_l = rho_h / r_rho;
 
         //kinimatic viscosity
         LBMReal nu_h = nuL;
         //LBMReal nu_l = nuG;
         //#dynamic viscosity
-        //LBMReal mu_h = rho_h * nu_h;
+        LBMReal mu_h = rho_h * nu_h;
         
         //gravity
         LBMReal g_y = Re * Re * nu_h * nu_h / (D*D*D);
@@ -174,7 +174,7 @@ void run(string configname)
         SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
         noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseNoSlipBCAlgorithm()));
         SPtr<BCAdapter> slipBCAdapter(new SlipBCAdapter());
-        slipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseSlipBCAlgorithm()));
+        noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseSlipBCAlgorithm()));
         //////////////////////////////////////////////////////////////////////////////////
         // BC visitor
         MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
@@ -197,55 +197,52 @@ void run(string configname)
         SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(grid, rSch, pathname, comm));
         //SPtr<MPIIOMigrationCoProcessor> rcp(new MPIIOMigrationCoProcessor(grid, rSch, metisVisitor, pathname, comm));
         //SPtr<MPIIOMigrationBECoProcessor> rcp(new MPIIOMigrationBECoProcessor(grid, rSch, metisVisitor, pathname, comm));
-        //rcp->setNu(nuLB);
-        //rcp->setNuLG(nuL, nuG);
+        // rcp->setNu(nuLB);
+       //  rcp->setNuLG(nuL, nuG);
         //rcp->setDensityRatio(densityRatio);
 
         rcp->setLBMKernel(kernel);
         rcp->setBCProcessor(bcProc);
         //////////////////////////////////////////////////////////////////////////
-        // bounding box
-        double g_minX1 = boundingBox[0];
-        double g_minX2 = boundingBox[2];
-        double g_minX3 = boundingBox[4];
-
-        double g_maxX1 = boundingBox[1];
-        double g_maxX2 = boundingBox[3];
-        double g_maxX3 = boundingBox[5];
-
-        double dx2 = 2.0 * dx;
-        GbCuboid3DPtr wallXmin(
-            new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_minX1, g_maxX2 + dx2, g_maxX3 + dx2));
-        GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
-        GbCuboid3DPtr wallXmax(
-            new GbCuboid3D(g_maxX1, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
-        GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
-
-        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> wallXminInt(new D3Q27Interactor(wallXmin, grid, slipBCAdapter, Interactor3D::SOLID));
-        SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, slipBCAdapter, Interactor3D::SOLID));
-
-        SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
-        SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
 
         if (newStart) {
 
+            // bounding box
+            double g_minX1 = boundingBox[0];
+            double g_minX2 = boundingBox[2];
+            double g_minX3 = boundingBox[4];
+
+            double g_maxX1 = boundingBox[1];
+            double g_maxX2 = boundingBox[3];
+            double g_maxX3 = boundingBox[5];
+
             // geometry
             SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
             if (myid == 0)
                 GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube",
-                                           WbWriterVtkXmlBinary::getInstance());
+                    WbWriterVtkXmlBinary::getInstance());
+
+
 
             GenBlocksGridVisitor genBlocks(gridCube);
             grid->accept(genBlocks);
 
+            double dx2 = 2.0 * dx;
+            GbCuboid3DPtr wallXmin(new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_minX1, g_maxX2 + dx2, g_maxX3 + dx2));
+            GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
+            GbCuboid3DPtr wallXmax(new GbCuboid3D(g_maxX1, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
+            GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
+
+            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> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+            SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
 
+            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));
@@ -352,28 +349,7 @@ void run(string configname)
             }
 
             rcp->restart((int)restartStep);
-            //grid->setTimeStep(restartStep);
-
-            //rcp->readBlocks((int)restartStep);
-            //rcp->readDataSet((int)restartStep);
-            //rcp->readBoundaryConds((int)restartStep);
-            //grid->setTimeStep((int)restartStep);
-
-            //SetBcBlocksBlockVisitor v2(wallXminInt);
-            //grid->accept(v2);
-            //wallXminInt->initInteractor();
-
-            //SetBcBlocksBlockVisitor v3(wallXmaxInt);
-            //grid->accept(v3);
-            //wallXmaxInt->initInteractor();
-
-            //SetBcBlocksBlockVisitor v4(wallYminInt);
-            //grid->accept(v4);
-            //wallYminInt->initInteractor();
-
-            //SetBcBlocksBlockVisitor v1(wallYmaxInt);
-            //grid->accept(v1);
-            //wallYmaxInt->initInteractor();
+            grid->setTimeStep(restartStep);
 
             if (myid == 0)
                 UBLOG(logINFO, "Restart - end");
@@ -416,9 +392,8 @@ void run(string configname)
 
         SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
             grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
-        //if(grid->getTimeStep() == 0) 
-            //pp->process(0);
-        pp->process(grid->getTimeStep());
+        if(grid->getTimeStep() == 0) 
+            pp->process(0);
 
         SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
         SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
index 0dc57aaba..df437b553 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
@@ -47,8 +47,8 @@ MultiphasePressureFilterLBMKernel::MultiphasePressureFilterLBMKernel() { this->c
 //////////////////////////////////////////////////////////////////////////
 void MultiphasePressureFilterLBMKernel::initDataSet()
 {
-	SPtr<DistributionArray3D> f(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9));
-	SPtr<DistributionArray3D> h(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9)); // For phase-field
+	SPtr<DistributionArray3D> f(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
+	SPtr<DistributionArray3D> h(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0)); // For phase-field
 
 	//SPtr<PhaseFieldArray3D> divU1(new PhaseFieldArray3D(            nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
 	CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure(new  CbArray3D<LBMReal, IndexerX3X2X1>(    nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
@@ -58,7 +58,7 @@ void MultiphasePressureFilterLBMKernel::initDataSet()
 	//dataSet->setPhaseField(divU1);
 	dataSet->setPressureField(pressure);
 
-	phaseField = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.0));
+	phaseField = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
 
 	divU = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
 }
@@ -156,6 +156,21 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 	nonLocalDistributionsH1 = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getNonLocalDistributions();
 	zeroDistributionsH1     = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getZeroDistributions();
 
+
+	/*std::vector<double> doubleValuesArrayFtemp; // double-values (arrays of f's) in all blocks  Fdistribution
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), localDistributionsF->getDataVector().begin(), localDistributionsF->getDataVector().end());
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), nonLocalDistributionsF->getDataVector().begin(), nonLocalDistributionsF->getDataVector().end());
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), zeroDistributionsF->getDataVector().begin(), zeroDistributionsF->getDataVector().end());
+
+	for (int i = 0; i < doubleValuesArrayFtemp.size(); i++)
+	{
+		if (UbMath::isNaN(doubleValuesArrayFtemp[i]) || UbMath::isInfinity(doubleValuesArrayFtemp[i]))
+			std::cout << "doubleValuesArrayFtemp[i] is NAN  " << i << std::endl;
+	}*/
+	//std::cout << "localDistributionsF=" << localDistributionsF->getNX1() << " " << localDistributionsF->getNX2() << " " << localDistributionsF->getNX3() << " " << localDistributionsF->getNX4() << std::endl;
+	//std::cout << "nonLocalDistributionsF=" << nonLocalDistributionsF->getNX1() << " " << nonLocalDistributionsF->getNX2() << " " << nonLocalDistributionsF->getNX3() << " " << nonLocalDistributionsF->getNX4() << std::endl;
+	//std::cout << "zeroDistributionsF=" << zeroDistributionsF->getNX1() << " " << zeroDistributionsF->getNX2() << " " << zeroDistributionsF->getNX3() << std::endl;
+
 	CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure = dataSet->getPressureField();
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
@@ -171,6 +186,12 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 	int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
 	int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
 
+//	std::cout << "MultiphasePressureFilterLBMKernel::calculate bcArrayMaxX1=" << bcArrayMaxX1 << " , bcArrayMaxX2=" << bcArrayMaxX2 << ", bcArrayMaxX3=" << bcArrayMaxX3 << std::endl;
+//	std::cout << "MultiphasePressureFilterLBMKernel::calculate minX1=" << minX1 << " , minX2=" << minX2 << ", minX3=" << minX3 << std::endl;
+//	std::cout << "MultiphasePressureFilterLBMKernel::calculate maxX1=" << maxX1 << " , maxX2=" << maxX2 << ", maxX3=" << maxX3 << std::endl;
+//	std::cout << "dataSetParamStr1=" << localDistributionsF->getNX1() << " " << localDistributionsF->getNX2() << " " << localDistributionsF->getNX3() << " " << localDistributionsF->getNX4() << std::endl;
+//	std::cout << "nx=" << dataSet->getFdistributions()->getNX1() << " " << dataSet->getFdistributions()->getNX2() << " " << dataSet->getFdistributions()->getNX3() << std::endl;
+
 	for (int x3 = minX3-ghostLayerWidth; x3 < maxX3+ghostLayerWidth; x3++) {
 		for (int x2 = minX2-ghostLayerWidth; x2 < maxX2+ghostLayerWidth; x2++) {
 			for (int x1 = minX1-ghostLayerWidth; x1 < maxX1+ghostLayerWidth; x1++) {
@@ -243,6 +264,34 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 
 					mfbbb = (*this->zeroDistributionsF)(x1, x2, x3);
 					
+					/*if (UbMath::isNaN(mfcbb) || UbMath::isInfinity(mfcbb)) std::cout << "mfcbb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbcb) || UbMath::isInfinity(mfbcb)) std::cout << "mfbcb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbbc) || UbMath::isInfinity(mfbbc)) std::cout << "mfbbc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfccb) || UbMath::isInfinity(mfccb)) std::cout << "mfccb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfacb) || UbMath::isInfinity(mfacb)) std::cout << "mfacb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcbc) || UbMath::isInfinity(mfcbc)) std::cout << "mfcbc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfabc) || UbMath::isInfinity(mfabc)) std::cout << "mfabc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbcc) || UbMath::isInfinity(mfbcc)) std::cout << "mfbcc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbac) || UbMath::isInfinity(mfbac)) std::cout << "mfbac is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfccc) || UbMath::isInfinity(mfccc)) std::cout << "mfccc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfacc) || UbMath::isInfinity(mfacc)) std::cout << "mfacc is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcac) || UbMath::isInfinity(mfcac)) std::cout << "mfcac is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfaac) || UbMath::isInfinity(mfaac)) std::cout << "mfaac is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfabb) || UbMath::isInfinity(mfabb)) std::cout << "mfabb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbab) || UbMath::isInfinity(mfbab)) std::cout << "mfbab is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbba) || UbMath::isInfinity(mfbba)) std::cout << "mfbba is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfaab) || UbMath::isInfinity(mfaab)) std::cout << "mfaab is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcab) || UbMath::isInfinity(mfcab)) std::cout << "mfcab is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfaba) || UbMath::isInfinity(mfaba)) std::cout << "mfaba is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcba) || UbMath::isInfinity(mfcba)) std::cout << "mfcba is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbaa) || UbMath::isInfinity(mfbaa)) std::cout << "mfbaa is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbca) || UbMath::isInfinity(mfbca)) std::cout << "mfbca is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfaaa) || UbMath::isInfinity(mfaaa)) std::cout << "mfaaa is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcaa) || UbMath::isInfinity(mfcaa)) std::cout << "mfcaa is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfaca) || UbMath::isInfinity(mfaca)) std::cout << "mfaca is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfcca) || UbMath::isInfinity(mfcca)) std::cout << "mfcca is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;
+					if (UbMath::isNaN(mfbbb) || UbMath::isInfinity(mfbbb)) std::cout << "mfbbb is NAN  " << x1 << " " << x2 << " " << x3 << std::endl;*/
+
 					LBMReal rhoH = 1.0;
 					LBMReal rhoL = 1.0 / densityRatio;
 
@@ -1200,7 +1249,8 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 					(*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac   ;//* rho * c1o3;
 
 					(*this->zeroDistributionsF)(x1, x2, x3) = mfbbb;// *rho* c1o3;
-																																		// !Old Kernel
+									
+																			// !Old Kernel
 /////////////////////  P H A S E - F I E L D   S O L V E R
 ////////////////////////////////////////////
 /////CUMULANT PHASE-FIELD
@@ -1462,6 +1512,23 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 			}
 		}
 	}
+
+	/*std::vector<double> doubleValuesArrayFtemp; // double-values (arrays of f's) in all blocks  Fdistribution
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), *this->localDistributionsF->getDataVector().begin(),
+		*this->localDistributionsF->getDataVector().end());
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), *this->nonLocalDistributionsF->getDataVector().begin(),
+		*this->nonLocalDistributionsF->getDataVector().end());
+	doubleValuesArrayFtemp.insert(doubleValuesArrayFtemp.end(), *this->zeroDistributionsF->getDataVector().begin(),
+		*this->zeroDistributionsF->getDataVector().end());
+
+	for (int i = 0; i < doubleValuesArrayFtemp.size(); i++)
+	{
+		if (UbMath::isNaN(doubleValuesArrayFtemp[i]) || UbMath::isInfinity(doubleValuesArrayFtemp[i]))
+			std::cout << "doubleValuesArrayFtemp[i] is NAN  " << i << std::endl;
+	}*/
+
+	//std::cout << "MPIIOMigrationCoProcessor::readDataSet posle NAN ni4ego net" << std::endl;*/
+
 }
 //////////////////////////////////////////////////////////////////////////
 
-- 
GitLab