Basic subgrid construction, coupling automatization

1 files changed, 83 insertions, 52 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index f551820..9f92c2a 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -197,6 +197,9 @@ void getResults(const std::string& prefix,
   }
 }
 
+template <typename T, template<typename> class DESCRIPTOR> class FineCoupler;
+template <typename T, template<typename> class DESCRIPTOR> class CoarseCoupler;
+
 template <typename T, template<typename> class DESCRIPTOR>
 class Grid {
 private:
@@ -206,6 +209,11 @@ private:
   std::unique_ptr<SuperGeometry2D<T>>           _geometry;
   std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice;
 
+  std::vector<std::unique_ptr<Grid<T,DESCRIPTOR>>> _fineGrids;
+
+  std::vector<std::unique_ptr<FineCoupler<T,DESCRIPTOR>>> _fineCouplers;
+  std::vector<std::unique_ptr<CoarseCoupler<T,DESCRIPTOR>>> _coarseCouplers;
+
 public:
   static std::unique_ptr<Grid<T,DESCRIPTOR>> make(IndicatorF2D<T>& domainF, int resolution, T tau, int re)
   {
@@ -274,15 +282,78 @@ public:
     return 1./getScalingFactor();
   }
 
-  std::unique_ptr<Grid<T,DESCRIPTOR>> refine(IndicatorF2D<T>& domainF)
+  void collideAndStream()
   {
-    return std::unique_ptr<Grid<T,DESCRIPTOR>>(
-             new Grid<T,DESCRIPTOR>(
-               domainF,
-               2*getConverter().getResolution(),
-               2.0*getConverter().getLatticeRelaxationTime() - 0.5,
-               getConverter().getReynoldsNumber()
-             ));
+    for ( auto& fineCoupler : _fineCouplers ) {
+      fineCoupler->store();
+    }
+
+    this->getSuperLattice().collideAndStream();
+
+    for ( auto& fineGrid : _fineGrids ) {
+      fineGrid->getSuperLattice().collideAndStream();
+    }
+
+    for ( auto& fineCoupler : _fineCouplers ) {
+      fineCoupler->interpolate();
+      fineCoupler->couple();
+    }
+
+    for ( auto& fineGrid : _fineGrids ) {
+      fineGrid->getSuperLattice().collideAndStream();
+    }
+
+    for ( auto& fineCoupler : _fineCouplers ) {
+      fineCoupler->store();
+      fineCoupler->couple();
+    }
+
+    for ( auto& coarseCoupler : _coarseCouplers ) {
+      coarseCoupler->couple();
+    }
+  }
+
+  Grid<T,DESCRIPTOR>* refine(Vector<T,2> origin, Vector<T,2> extend)
+  {
+    IndicatorCuboid2D<T> fineCuboid(extend, origin);
+    _fineGrids.emplace_back(
+      new Grid<T,DESCRIPTOR>(
+        fineCuboid,
+        2*getConverter().getResolution(),
+        2.0*getConverter().getLatticeRelaxationTime() - 0.5,
+        getConverter().getReynoldsNumber()
+      ));
+    Grid<T,DESCRIPTOR>* const fineGrid = _fineGrids.back().get();
+
+    const T coarseDeltaX = getConverter().getPhysDeltaX();
+
+    _fineCouplers.emplace_back(
+      new FineCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin, Vector<T,2> {0,extend[1]}));
+    _fineCouplers.emplace_back(
+      new FineCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + Vector<T,2> {extend[0],0}, Vector<T,2> {0,extend[1]}));
+    _fineCouplers.emplace_back(
+      new FineCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + Vector<T,2>(0,extend[1]), Vector<T,2> {extend[0],0}));
+    _fineCouplers.emplace_back(
+      new FineCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin, Vector<T,2> {extend[0],0}));
+
+    _coarseCouplers.emplace_back(
+      new CoarseCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
+    _coarseCouplers.emplace_back(
+      new CoarseCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + Vector<T,2>(extend[0]-coarseDeltaX,coarseDeltaX), Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
+    _coarseCouplers.emplace_back(
+      new CoarseCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + Vector<T,2> {coarseDeltaX,extend[1]-coarseDeltaX}, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));
+    _coarseCouplers.emplace_back(
+      new CoarseCoupler<T,DESCRIPTOR>(
+        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));
+
+    return fineGrid;
   }
 };
 
@@ -619,12 +690,11 @@ int main(int argc, char* argv[])
   Vector<T,2> coarseExtend { lx, ly  };
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
 
-  Vector<T,2> fineOrigin   { 0.25 * lx, 0.3 };
-  Vector<T,2> fineExtend   { 0.5  * lx, 0.4 };
-  IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
+  Vector<T,2> fineOrigin   { 0.25 * lx, 0.2 };
+  Vector<T,2> fineExtend   { 0.6  * lx, 0.6 };
 
   auto coarseGrid = Grid<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
-  auto fineGrid   = coarseGrid->refine(fineCuboid);
+  auto fineGrid   = coarseGrid->refine(fineOrigin, fineExtend);
 
   prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
   prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry());
@@ -675,52 +745,13 @@ int main(int argc, char* argv[])
     residuum);
   timer.start();
 
-  FineCoupler<T,DESCRIPTOR> fineCoupler0(*coarseGrid, *fineGrid, fineOrigin, Vector<T,2> {0,fineExtend[1]});
-  FineCoupler<T,DESCRIPTOR> fineCoupler1(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2> {fineExtend[0],0}, Vector<T,2> {0,fineExtend[1]});
-  FineCoupler<T,DESCRIPTOR> fineCoupler2(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2>(0,fineExtend[1]), Vector<T,2> {fineExtend[0],0});
-  FineCoupler<T,DESCRIPTOR> fineCoupler3(*coarseGrid, *fineGrid, fineOrigin, Vector<T,2> {fineExtend[0],0});
-
-  CoarseCoupler<T,DESCRIPTOR> coarseCoupler0(*coarseGrid, *fineGrid, fineOrigin + coarseDeltaX, Vector<T,2> {0,fineExtend[1]-2*coarseDeltaX});
-  CoarseCoupler<T,DESCRIPTOR> coarseCoupler1(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2>(fineExtend[0]-coarseDeltaX,coarseDeltaX), Vector<T,2> {0,fineExtend[1]-2*coarseDeltaX});
-  CoarseCoupler<T,DESCRIPTOR> coarseCoupler2(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2> {coarseDeltaX,fineExtend[1]-coarseDeltaX}, Vector<T,2> {fineExtend[0]-2*coarseDeltaX,0});
-  CoarseCoupler<T,DESCRIPTOR> coarseCoupler3(*coarseGrid, *fineGrid, fineOrigin + coarseDeltaX, Vector<T,2> {fineExtend[0]-2*coarseDeltaX,0});
-
   for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
       clout << "Simulation converged." << endl;
       break;
     }
 
-    fineCoupler0.store();
-    fineCoupler1.store();
-    fineCoupler2.store();
-    fineCoupler3.store();
-    coarseGrid->getSuperLattice().collideAndStream();
-
-    fineGrid->getSuperLattice().collideAndStream();
-    fineCoupler0.interpolate();
-    fineCoupler0.couple();
-    fineCoupler1.interpolate();
-    fineCoupler1.couple();
-    fineCoupler2.interpolate();
-    fineCoupler2.couple();
-    fineCoupler3.interpolate();
-    fineCoupler3.couple();
-
-    fineGrid->getSuperLattice().collideAndStream();
-    fineCoupler0.store();
-    fineCoupler0.couple();
-    fineCoupler1.store();
-    fineCoupler1.couple();
-    fineCoupler2.store();
-    fineCoupler2.couple();
-    fineCoupler3.store();
-    fineCoupler3.couple();
-
-    coarseCoupler0.couple();
-    coarseCoupler1.couple();
-    coarseCoupler2.couple();
-    coarseCoupler3.couple();
+    coarseGrid->collideAndStream();
 
     getResults(
       "coarse_",