From c7955c90d820a64fb6315d35ff510660fcfea89c Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Sun, 30 Dec 2018 22:07:47 +0100
Subject: Test order 4 interpolation of full distributions

i.e. not only the moments and non-equilibrium parts.
---
 apps/adrian/poiseuille2d/poiseuille2d.cpp | 45 +++++++++++++++++++++++++------
 1 file changed, 37 insertions(+), 8 deletions(-)

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 43f8631..583282e 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -48,7 +48,6 @@ const T Re = 10.;             // Reynolds number
 const T maxPhysT = 40.;       // max. simulation time in s, SI unit
 const T physInterval = 0.25;  // interval for the convergence check in s
 const T residuum = 1e-5;      // residuum for the convergence check
-const T tuner = 0.99;         // for partialSlip only: 0->bounceBack, 1->freeSlip
 
 
 void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
@@ -59,7 +58,7 @@ void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
   clout << "Prepare Geometry ..." << std::endl;
 
   superGeometry.rename(0,2);
-
+  // Set material number for bounce back boundaries
   superGeometry.rename(2,1,0,1);
 
   Vector<T,2> extend;
@@ -193,7 +192,7 @@ void getResults(const std::string& prefix,
     vtmWriter.createMasterFile();
   }
 
-  if (iT%50==0) {
+  if (iT%100==0) {
     vtmWriter.write(iT);
   }
 
@@ -218,7 +217,7 @@ public:
   Grid(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
     _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
                  resolution, // resolution: number of voxels per charPhysL
-                 tau,        // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
+                 tau,        // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
                  T{1},       // charPhysLength: reference length of simulation geometry
                  T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
                  T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
@@ -265,6 +264,16 @@ public:
   }
 };
 
+T order4interpolation(T fm3, T fm1, T f1, T f3)
+{
+  return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
+}
+
+T order3interpolation(T fm1, T f1, T f3)
+{
+  return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
+}
+
 template <typename T, template<typename> class DESCRIPTOR>
 void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
 {
@@ -272,16 +281,36 @@ void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
   auto& fineLattice   = fine.getSuperLattice().getBlockLattice(0);
 
   const int x = coarseLattice.getNx()-2;
+
   for (int y=0; y < coarseLattice.getNy(); ++y) {
     for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
       fineLattice.get(0, 2*y)[i] = coarseLattice.get(x,y)[i];
     }
   }
-  for (int y=0; y < coarseLattice.getNy()-1; ++y) {
+
+  for (int y=2; y < coarseLattice.getNy()-2; ++y) {
     for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      fineLattice.get(0,1+2*y)[i] = 0.5 * (coarseLattice.get(x,y)[i] + coarseLattice.get(x,y+1)[i]);
+      fineLattice.get(0,1+2*y)[i] = order4interpolation(
+        coarseLattice.get(x,y-1)[i],
+        coarseLattice.get(x,y)[i],
+        coarseLattice.get(x,y+1)[i],
+        coarseLattice.get(x,y+2)[i]
+      );
     }
   }
+
+  for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
+    fineLattice.get(0,1)[i] = order3interpolation(
+      coarseLattice.get(x,0)[i],
+      coarseLattice.get(x,1)[i],
+      coarseLattice.get(x,2)[i]
+    );
+    fineLattice.get(0,fineLattice.getNy()-2)[i] = order3interpolation(
+      coarseLattice.get(x,coarseLattice.getNy()-1)[i],
+      coarseLattice.get(x,coarseLattice.getNy()-2)[i],
+      coarseLattice.get(x,coarseLattice.getNy()-3)[i]
+    );
+  }
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -362,7 +391,7 @@ int main(int argc, char* argv[])
     coarseGrid.getConverter().getLatticeTime(maxPhysT),
     coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
   util::ValueTracer<T> converge(
-    coarseGrid.getConverter().getLatticeTime(physInterval),
+    fineGrid.getConverter().getLatticeTime(physInterval),
     residuum);
   timer.start();
 
@@ -396,7 +425,7 @@ int main(int argc, char* argv[])
       timer,
       converge.hasConverged());
 
-    converge.takeValue(coarseGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
+    converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
   }
 
   timer.stop();
-- 
cgit v1.2.3