From b0cc9ac8f20b09b86ddef06b4188982e82ab8a57 Mon Sep 17 00:00:00 2001
From: "AMATERASU\\geier" <geier@irmb.tu-bs.de>
Date: Fri, 18 Jun 2021 15:52:21 +0200
Subject: [PATCH] output fixed for prototype velocity formulation without
 viscosity correction

---
 .../WriteMultiphaseQuantitiesCoProcessor.cpp  | 58 +++++++++++++------
 ...woPhaseFieldsVelocityCumulantLBMKernel.cpp |  2 +-
 ...eTwoPhaseFieldsVelocityCumulantLBMKernel.h |  7 ++-
 3 files changed, 47 insertions(+), 20 deletions(-)

diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
index 70d2f2b6c..05ec0c36c 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -36,6 +36,7 @@
 #include "LBMKernel.h"
 #include <string>
 #include <vector>
+#include "MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.h"
 
 #include "BCArray3D.h"
 #include "Block3D.h"
@@ -145,7 +146,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 {
     using namespace D3Q27System;
     using namespace UbMath;
-
+    SPtr<LBMKernel> kernel = dynamicPointerCast<LBMKernel>(block->getKernel());
     //double level   = (double)block->getLevel();
 
     // Diese Daten werden geschrieben:
@@ -156,10 +157,11 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
     datanames.push_back("Vz");
     datanames.push_back("P1");
     datanames.push_back("Phi2");
+    if (dynamicPointerCast<MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel>(kernel)->pressure) datanames.push_back("Pressure");
 
     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();
@@ -227,6 +229,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                             ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
                 }
                     else { (*phaseField2)(ix1, ix2, ix3) = 999.0; }
+                    
                 }
             }
         }
@@ -344,23 +347,40 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     // 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);
+                    if (dynamicPointerCast<MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel>(kernel)->pressure) {
+                        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])) ;
 
-                    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);
+                        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])) ;
 
-                    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);
+                        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]));
+
+                    }
+                    else {
+                        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])) +
@@ -409,6 +429,10 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     data[index++].push_back(vx3);
                     data[index++].push_back(p1);
                     data[index++].push_back(phi2[REST]);
+                    if (dynamicPointerCast<MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel>(kernel)->pressure) {
+                        data[index++].push_back((*dynamicPointerCast<MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel>(kernel)->pressure)(ix1, ix2, ix3));
+                    }
+                   // else { data[index++].push_back(999); }
                 }
             }
         }
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
index fd0b7f2a1..db01b3edc 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
@@ -52,7 +52,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::initDataSet()
     SPtr<DistributionArray3D> h2(new D3Q27EsoTwist3DSplittedVector(nx[0] + 2, nx[1] + 2, nx[2] + 2, -999.9)); // For phase-field
     SPtr<PhaseFieldArray3D> divU(new PhaseFieldArray3D(nx[0] + 2, nx[1] + 2, nx[2] + 2, 0.0));
 	 pressure= CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new  CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 2, nx[1] + 2, nx[2] + 2, 0.0));
-	 pressureOld = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new  CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 2, nx[1] + 2, nx[2] + 2, 0.0));
+	// pressureOld = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new  CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 2, nx[1] + 2, nx[2] + 2, 0.0));
     dataSet->setFdistributions(f);
     dataSet->setHdistributions(h); // For phase-field
     dataSet->setH2distributions(h2); // For phase-field
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.h
index 67232529c..7e14cc3c3 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.h
@@ -56,6 +56,10 @@ public:
    void forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
    void backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
 
+   ///refactor
+   CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure;
+   //CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressureOld;
+
    double getCalculationTime() override { return .0; }
 protected:
    virtual void initDataSet();
@@ -75,8 +79,7 @@ protected:
    CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr zeroDistributionsH2;
 
    //CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr   phaseField;
-   CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure;
-   CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressureOld;
+
 
    LBMReal h  [D3Q27System::ENDF+1];
    LBMReal h2[D3Q27System::ENDF + 1];
-- 
GitLab