From 8a2992826d153618765e8c406962c3dddd6cdbcd Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Tue, 29 Jan 2019 16:37:39 +0100
Subject: [PATCH] add pressure correction with Laplace operator

---
 source/Applications/ConvectionOfVortex/cov.cpp | 17 +++++++++--------
 source/Applications/VirtualFluids.h            |  1 +
 .../CompressibleOffsetInterpolationProcessor.h |  1 -
 ...ssibleOffsetMomentsInterpolationProcessor.h |  1 -
 ...setSquarePressureInterpolationProcessor.cpp | 18 ++++++++++++++----
 ...ffsetSquarePressureInterpolationProcessor.h |  1 -
 6 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/source/Applications/ConvectionOfVortex/cov.cpp b/source/Applications/ConvectionOfVortex/cov.cpp
index 37e66708c..a98b48f17 100644
--- a/source/Applications/ConvectionOfVortex/cov.cpp
+++ b/source/Applications/ConvectionOfVortex/cov.cpp
@@ -21,7 +21,7 @@ void run()
       //////////////////////////////////////////////////////////////////////////
       //DLR-F16 test
       //dx_coarse = 0.003 mm
-      //string  pathname = "d:/temp/ConvectionOfVortex_0.003";
+      //string  pathname = "d:/temp/ConvectionOfVortex_0.003_square";
       //int     endTime = 20;
       //double  outTime = 10;
       //LBMReal dx =  0.003;
@@ -38,8 +38,8 @@ void run()
       ////////////////////////////////////////////////////////////////////////////
       //dx_coarse = 0.00075 mm
       string  pathname = "d:/temp/ConvectionOfVortex_0.00075_moments";
-      double  endTime = 80;
-      double  outTime = 80;
+      double  endTime = 160;
+      double  outTime = 160;
       LBMReal dx =  0.00075;
       LBMReal rhoLB = 0.0;
       LBMReal nuLB = 8.66025e-6*4.0;
@@ -115,6 +115,7 @@ void run()
       //SPtr<InterpolationProcessor> iProcessor(new CompressibleOffsetInterpolationProcessor());
       SPtr<InterpolationProcessor> iProcessor(new CompressibleOffsetMomentsInterpolationProcessor());
       //dynamicPointerCast<CompressibleOffsetMomentsInterpolationProcessor>(iProcessor)->setBulkOmegaToOmega(true);
+      //SPtr<InterpolationProcessor> iProcessor(new CompressibleOffsetSquarePressureInterpolationProcessor());
       SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
 
       UBLOG(logINFO, "SetConnectorsBlockVisitor:start");
@@ -226,16 +227,16 @@ void run()
       SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
       std::shared_ptr<NUPSCounterCoProcessor> nupsCoProcessor(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
 
-      SPtr<UbScheduler> tavSch(new UbScheduler(1, 0, endTime));
-      SPtr<TimeAveragedValuesCoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, comm,
-         TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations));
-      tav->setWithGhostLayer(true);
+      //SPtr<UbScheduler> tavSch(new UbScheduler(1, 0, endTime));
+      //SPtr<TimeAveragedValuesCoProcessor> tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, comm,
+      //   TimeAveragedValuesCoProcessor::Density | TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations));
+      //tav->setWithGhostLayer(true);
 
       SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
       SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
       calculator->addCoProcessor(nupsCoProcessor);
       calculator->addCoProcessor(writeMQCoProcessor);
-      calculator->addCoProcessor(tav);
+      //calculator->addCoProcessor(tav);
 
       //omp_set_num_threads(1);
 
diff --git a/source/Applications/VirtualFluids.h b/source/Applications/VirtualFluids.h
index 95a3197ee..580b342a6 100644
--- a/source/Applications/VirtualFluids.h
+++ b/source/Applications/VirtualFluids.h
@@ -174,6 +174,7 @@
 #include <LBM/IncompressibleOffsetInterpolationProcessor.h>
 #include <LBM/CompressibleOffsetInterpolationProcessor.h>
 #include <LBM/CompressibleOffsetMomentsInterpolationProcessor.h>
+#include <LBM/CompressibleOffsetSquarePressureInterpolationProcessor.h>
 #include <LBM/InterpolationHelper.h>
 #include <LBM/InterpolationProcessor.h>
 //#include <LBM/D3Q27OffsetInterpolationProcessor.h>
diff --git a/source/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.h b/source/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.h
index d0dd40c30..2f3da5f03 100644
--- a/source/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.h
+++ b/source/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.h
@@ -10,7 +10,6 @@
 //////////////////////////////////////////////////////////////////////////
 
 class CompressibleOffsetInterpolationProcessor;
-typedef SPtr<CompressibleOffsetInterpolationProcessor> CompressibleOffsetInterpolationProcessorPtr;
 
 class CompressibleOffsetInterpolationProcessor : public InterpolationProcessor
 {
diff --git a/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h b/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
index b27e78663..f9876cf54 100644
--- a/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
+++ b/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
@@ -10,7 +10,6 @@
 //////////////////////////////////////////////////////////////////////////
 
 class CompressibleOffsetMomentsInterpolationProcessor;
-typedef SPtr<CompressibleOffsetMomentsInterpolationProcessor> CompressibleOffsetMomentsInterpolationProcessorPtr;
 
 class CompressibleOffsetMomentsInterpolationProcessor : public InterpolationProcessor
 {
diff --git a/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp b/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
index f272b5676..44be4234f 100644
--- a/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
+++ b/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
@@ -485,6 +485,11 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    LBMReal oP = OxxPyyPzzF;
 
    LBMReal rho  = press ;//+ (2.*axx*x+axy*y+axz*z+axyz*y*z+ax + 2.*byy*y+bxy*x+byz*z+bxyz*x*z+by + 2.*czz*z+cxz*x+cyz*y+cxyz*x*y+cz)/3.;
+
+   LBMReal laplaceRho = (xoff!=0.0 || yoff!=0.0 || zoff!= 0.0) ? 0.0 :(-3.0*(by*by+ax*ax+cz*cz)-6.0*(ay*bx+bz*cy+az*cx))*(1.0+rho);
+
+   rho=rho+laplaceRho*(3.0/16.0);
+
    LBMReal vx1  = a0 + 0.25*( xs*ax + ys*ay + zs*az) + 0.0625*(axx + xs*ys*axy + xs*zs*axz + ayy + ys*zs*ayz + azz) + 0.015625*(xs*ys*zs*axyz);
    LBMReal vx2  = b0 + 0.25*( xs*bx + ys*by + zs*bz) + 0.0625*(bxx + xs*ys*bxy + xs*zs*bxz + byy + ys*zs*byz + bzz) + 0.015625*(xs*ys*zs*bxyz);
    LBMReal vx3  = c0 + 0.25*( xs*cx + ys*cy + zs*cz) + 0.0625*(cxx + xs*ys*cxy + xs*zs*cxz + cyy + ys*zs*cyz + czz) + 0.015625*(xs*ys*zs*cxyz);
@@ -517,7 +522,7 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    LBMReal mfcaa = zeroReal;
    LBMReal mfaca = zeroReal;
 
-   mfaaa = press; // if drho is interpolated directly
+   mfaaa = rho; // if drho is interpolated directly
 
    LBMReal vx1Sq = vx1*vx1;
    LBMReal vx2Sq = vx2*vx2;
@@ -926,10 +931,15 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    LBMReal vx1  = a0;
    LBMReal vx2  = b0;
    LBMReal vx3  = c0;
-
+  
+   
    LBMReal rho = press ;//+ (ax+by+cz)/3.;
 
-   LBMReal eps_new = 2.;
+   LBMReal laplaceRho = (xoff!=0.0 || yoff!=0.0 || zoff!= 0.0) ? 0.0 :(-3.0*(by*by+ax*ax+cz*cz)-6.0*(ay*bx+bz*cy+az*cx))*(1.0+rho);
+
+   rho=rho-laplaceRho*0.25;
+
+   LBMReal eps_new = 2.0;
    LBMReal o  = omega;
    //bulk viscosity
    LBMReal oP = OxxPyyPzzC;
@@ -962,7 +972,7 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    LBMReal mfcaa = zeroReal;
    LBMReal mfaca = zeroReal;
 
-   mfaaa = press; // if drho is interpolated directly
+   mfaaa = rho; // if drho is interpolated directly
 
    LBMReal vx1Sq = vx1*vx1;
    LBMReal vx2Sq = vx2*vx2;
diff --git a/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.h b/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.h
index 21f88accb..48aed6665 100644
--- a/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.h
+++ b/source/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.h
@@ -10,7 +10,6 @@
 //////////////////////////////////////////////////////////////////////////
 
 class CompressibleOffsetSquarePressureInterpolationProcessor;
-typedef SPtr<CompressibleOffsetSquarePressureInterpolationProcessor> CompressibleOffsetMomentsInterpolationProcessorPtr;
 
 class CompressibleOffsetSquarePressureInterpolationProcessor : public InterpolationProcessor
 {
-- 
GitLab