Test basic coarse subset refinement

1 files changed, 28 insertions, 18 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 632fba5..6c50339 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -97,8 +97,10 @@ void prepareCoarseLattice(UnitConverter<T,DESCRIPTOR> const& converter,
   sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
   sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
   sLattice.defineDynamics(superGeometry, 3, &bulkDynamics); // inflow
+  sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); // outflow
 
   sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);
+  sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);
 
   const T Lx = converter.getLatticeLength(lx);
   const T Ly = converter.getLatticeLength(ly);
@@ -118,6 +120,8 @@ void prepareCoarseLattice(UnitConverter<T,DESCRIPTOR> const& converter,
 
   sLattice.defineRhoU(superGeometry, 3, rho, u);
   sLattice.iniEquilibrium(superGeometry, 3, rho, u);
+  sLattice.defineRhoU(superGeometry, 4, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 4, rho, u);
 
   sLattice.initialize();
 
@@ -138,9 +142,6 @@ void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
   sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
   sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
   sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
-  sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); // outflow
-
-  sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);
 
   const T Lx = converter.getLatticeLength(lx);
   const T Ly = converter.getLatticeLength(ly);
@@ -149,7 +150,7 @@ void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
 
   const T maxVelocity = converter.getCharLatticeVelocity();
   const T radius = ly/2;
-  std::vector<T> axisPoint{lx, ly/2};
+  std::vector<T> axisPoint{0, ly/2};
   std::vector<T> axisDirection{1, 0};
   Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);
 
@@ -158,9 +159,6 @@ void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
   sLattice.defineRhoU(superGeometry, 2, rho, u);
   sLattice.iniEquilibrium(superGeometry, 2, rho, u);
 
-  sLattice.defineRhoU(superGeometry, 4, rho, u);
-  sLattice.iniEquilibrium(superGeometry, 4, rho, u);
-
   sLattice.initialize();
 
   clout << "Prepare fine lattice ... OK" << std::endl;
@@ -610,11 +608,11 @@ int main(int argc, char* argv[])
   const T coarseDeltaX = 1./N;
 
   Vector<T,2> coarseOrigin { 0.0, 0.0 };
-  Vector<T,2> coarseExtend { 0.5 * lx, ly  };
+  Vector<T,2> coarseExtend { lx, ly  };
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
 
-  Vector<T,2> fineOrigin   { 0.5 * lx - coarseDeltaX, 0.0 };
-  Vector<T,2> fineExtend   { 0.5 * lx + coarseDeltaX, ly  };
+  Vector<T,2> fineOrigin   { 0.25 * lx, 0.0 };
+  Vector<T,2> fineExtend   { 0.5  * lx, ly };
   IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
 
   auto coarseGrid = Grid<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
@@ -623,6 +621,10 @@ int main(int argc, char* argv[])
   prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
   prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry());
 
+  IndicatorCuboid2D<T> overlapCuboid(Vector<T,2>{fineExtend[0]-2*coarseDeltaX, fineExtend[1]}, Vector<T,2>{fineOrigin[0]+coarseDeltaX, fineOrigin[1]});
+  coarseGrid->getSuperGeometry().rename(1,0,overlapCuboid);
+  coarseGrid->getSuperGeometry().rename(2,0,overlapCuboid);
+
   Dynamics<T, DESCRIPTOR>* coarseBulkDynamics;
   coarseBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
     coarseGrid->getConverter().getLatticeRelaxationFrequency(),
@@ -664,11 +666,13 @@ int main(int argc, char* argv[])
 
   const Vector<T,2> originC2F = fineOrigin;
   const Vector<T,2> extendC2F { 0, fineExtend[1] };
-  FineCoupler<T,DESCRIPTOR> fineCoupler(*coarseGrid, *fineGrid, originC2F, extendC2F);
+  FineCoupler<T,DESCRIPTOR> fineCoupler0(*coarseGrid, *fineGrid, fineOrigin, extendC2F);
+  FineCoupler<T,DESCRIPTOR> fineCoupler1(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2>(fineExtend[0],0), extendC2F);
 
   const Vector<T,2> originF2C { fineOrigin[0] + coarseDeltaX, fineOrigin[1] };
   const Vector<T,2> extendF2C { 0, fineExtend[1] };
-  CoarseCoupler<T,DESCRIPTOR> coarseCoupler(*coarseGrid, *fineGrid, originF2C, extendF2C);
+  CoarseCoupler<T,DESCRIPTOR> coarseCoupler0(*coarseGrid, *fineGrid, originF2C, extendF2C);
+  CoarseCoupler<T,DESCRIPTOR> coarseCoupler1(*coarseGrid, *fineGrid, fineOrigin + Vector<T,2>(fineExtend[0]-coarseDeltaX,0), extendF2C);
 
   for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
@@ -676,18 +680,24 @@ int main(int argc, char* argv[])
       break;
     }
 
-    fineCoupler.store();
+    fineCoupler0.store();
+    fineCoupler1.store();
     coarseGrid->getSuperLattice().collideAndStream();
 
     fineGrid->getSuperLattice().collideAndStream();
-    fineCoupler.interpolate();
-    fineCoupler.couple();
+    fineCoupler0.interpolate();
+    fineCoupler0.couple();
+    fineCoupler1.interpolate();
+    fineCoupler1.couple();
 
     fineGrid->getSuperLattice().collideAndStream();
-    fineCoupler.store();
-    fineCoupler.couple();
+    fineCoupler0.store();
+    fineCoupler0.couple();
+    fineCoupler1.store();
+    fineCoupler1.couple();
 
-    coarseCoupler.couple();
+    coarseCoupler0.couple();
+    coarseCoupler1.couple();
 
     getResults(
       "coarse_",