Bidirectional coupling without regard for correctness

Looks surprisingly well from a purely visual perspective.

1 files changed, 45 insertions, 27 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 9415359..43f8631 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -43,9 +43,9 @@ typedef double T;
 
 const T lx  = 2.;             // length of the channel
 const T ly  = 1.;             // height of the channel
-int N = 20;                   // resolution of the model
+const int N = 20;             // resolution of the model
 const T Re = 10.;             // Reynolds number
-const T maxPhysT = 20.;       // max. simulation time in s, SI unit
+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
@@ -187,14 +187,13 @@ void getResults(const std::string& prefix,
   vtmWriter.addFunctor(velocity);
   vtmWriter.addFunctor(pressure);
 
-  const int vtmIter  = converter.getLatticeTime(maxPhysT/100.);
   const int statIter = converter.getLatticeTime(maxPhysT/10.);
 
-  if ( iT==0 ) {
+  if (iT==0) {
     vtmWriter.createMasterFile();
   }
 
-  if (iT%vtmIter==0 || hasConverged) {
+  if (iT%50==0) {
     vtmWriter.write(iT);
   }
 
@@ -266,6 +265,39 @@ public:
   }
 };
 
+template <typename T, template<typename> class DESCRIPTOR>
+void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
+{
+  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
+  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 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]);
+    }
+  }
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+void coupleF2C(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
+{
+  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
+  auto& fineLattice   = fine.getSuperLattice().getBlockLattice(0);
+
+  const int x = coarseLattice.getNx()-1;
+  for (int y=0; y < coarseLattice.getNy(); ++y) {
+    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
+      coarseLattice.get(x,y)[i] = fineLattice.get(2,2*y)[i];
+    }
+  }
+}
+
 int main(int argc, char* argv[])
 {
   olbInit(&argc, &argv);
@@ -293,7 +325,6 @@ int main(int argc, char* argv[])
     Re);
 
   prepareGeometry(coarseGrid.getConverter(), coarseGrid.getSuperGeometry());
-
   prepareGeometry(fineGrid.getConverter(), fineGrid.getSuperGeometry());
 
   Dynamics<T, DESCRIPTOR>* coarseBulkDynamics;
@@ -335,31 +366,27 @@ int main(int argc, char* argv[])
     residuum);
   timer.start();
 
-  for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); iT += 2) {
+  for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
       clout << "Simulation converged." << endl;
-      getResults(
-        "coarse_",
-        coarseGrid.getSuperLattice(),
-        coarseGrid.getConverter(),
-        iT,
-        coarseGrid.getSuperGeometry(),
-        timer,
-        converge.hasConverged());
       break;
     }
 
     coarseGrid.getSuperLattice().collideAndStream();
+    fineGrid.getSuperLattice().collideAndStream();
+    coupleC2F(coarseGrid, fineGrid);
+    fineGrid.getSuperLattice().collideAndStream();
+    coupleC2F(coarseGrid, fineGrid);
+    coupleF2C(coarseGrid, fineGrid);
+
     getResults(
       "coarse_",
       coarseGrid.getSuperLattice(),
       coarseGrid.getConverter(),
-      iT / 2,
+      iT,
       coarseGrid.getSuperGeometry(),
       timer,
       converge.hasConverged());
-
-    fineGrid.getSuperLattice().collideAndStream();
     getResults(
       "fine_",
       fineGrid.getSuperLattice(),
@@ -368,15 +395,6 @@ int main(int argc, char* argv[])
       fineGrid.getSuperGeometry(),
       timer,
       converge.hasConverged());
-    fineGrid.getSuperLattice().collideAndStream();
-    getResults(
-      "fine_",
-      fineGrid.getSuperLattice(),
-      fineGrid.getConverter(),
-      iT+1,
-      fineGrid.getSuperGeometry(),
-      timer,
-      converge.hasConverged());
 
     converge.takeValue(coarseGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
   }