Refine cylinder2d to stay stable even at very low resolutions

Somewhat fiddled together but works as expected. Refinement areas are
informed by Knudsen metric results. The key is to refine the outflow.
Sadly this requires refinement-overlap-crossing boundary conditions.

Interface for detailed (i.e. cell-relative) grid refinement control
still has room for improvement. Mixing physical indicators and cell
alignment requirements is problematic.

1 files changed, 49 insertions, 34 deletions

diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp
index ffab6a0..bf25417 100644
--- a/apps/adrian/cylinder2d/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/cylinder2d.cpp
@@ -39,11 +39,11 @@ typedef double T;
 /// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
 const T cylinderD = 1.0;
 
-const int N = 20;              // resolution of the model
+const int N = 5;               // resolution of the model
 const T lx  = 22  * cylinderD; // length of the channel
 const T ly  = 4.1 * cylinderD; // height of the channel
 const T Re  = 100.;            // Reynolds number
-const T uLat = 0.1;            // lattice velocity
+const T uLat = 0.05;           // lattice velocity
 const T maxPhysT = 60.;        // max. simulation time in s, SI unit
 
 
@@ -67,13 +67,13 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
     IndicatorCuboid2D<T> lowerWall(wallExtend, wallOrigin);
     sGeometry.rename(1,2,lowerWall);
 
-    IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin + Vector<T,2> {0,ly});
+    IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin + Vector<T,2> {0,ly-physSpacing/2});
     sGeometry.rename(1,2,upperWall);
   }
 
   // Set material number for inflow and outflow
   {
-    const Vector<T,2> extend { physSpacing/2,  ly};
+    const Vector<T,2> extend { physSpacing/2,  ly-physSpacing/2};
     const Vector<T,2> origin {-physSpacing/4, -physSpacing/4};
 
     IndicatorCuboid2D<T> inflow(extend, origin);
@@ -216,27 +216,58 @@ int main(int argc, char* argv[])
     coarseGrid.getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseBoundaryCondition(coarseGrid.getSuperLattice());
+  createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseBoundaryCondition);
 
-  const Vector<T,2> fineExtend {16*cylinderD, 3.75*cylinderD};
-  const Vector<T,2> fineOrigin {0.75*cylinderD, (ly-fineExtend[1])/2};
+  const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
+
+  const Vector<T,2> fineExtend {8*cylinderD, ly-2*coarseDeltaX};
+  const Vector<T,2> fineOrigin {0.75*cylinderD, coarseDeltaX};
 
   auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
   prepareGeometry(fineGrid);
 
+  const Vector<T,2> fineExtendB {10*coarseDeltaX, ly};
+  const Vector<T,2> fineOriginB {lx-10*coarseDeltaX, 0};
+
+  auto& fineOutflowGrid = coarseGrid.refine(fineOriginB, fineExtendB, false);
+  prepareGeometry(fineOutflowGrid);
+
+  {
+    const Vector<T,2> origin = fineOutflowGrid.getOrigin();
+    const Vector<T,2> extend = fineOutflowGrid.getExtend();
+
+    const Vector<T,2> extendY = {0,extend[1]};
+
+    coarseGrid.addFineCoupling(fineOutflowGrid, origin, extendY);
+    coarseGrid.addCoarseCoupling(fineOutflowGrid, origin + Vector<T,2> {coarseDeltaX,0}, extendY);
+
+    IndicatorCuboid2D<T> refined(extend, origin + Vector<T,2> {2*coarseDeltaX,0});
+
+    coarseGrid.getSuperGeometry().rename(1,0,refined);
+    coarseGrid.getSuperGeometry().rename(2,0,refined);
+    coarseGrid.getSuperGeometry().rename(4,0,refined);
+  }
+
   disableRefinedArea(coarseGrid, fineGrid);
 
   BGKdynamics<T, DESCRIPTOR> fineBulkDynamics(
     fineGrid.getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineBoundaryCondition(fineGrid.getSuperLattice());
+  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineBoundaryCondition);
 
-  prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition);
+  BGKdynamics<T, DESCRIPTOR> fineOutflowBulkDynamics(
+    fineOutflowGrid.getConverter().getLatticeRelaxationFrequency(),
+    instances::getBulkMomenta<T, DESCRIPTOR>());
+
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineOutflowBoundaryCondition(fineOutflowGrid.getSuperLattice());
+  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineOutflowBoundaryCondition);
+
+  prepareLattice(coarseGrid, coarseBulkDynamics, coarseBoundaryCondition);
 
-  const Vector<T,2> fineExtend2 {11*cylinderD, 3*cylinderD};
+  const Vector<T,2> fineExtend2 {4*cylinderD, 2*cylinderD};
   const Vector<T,2> fineOrigin2 {1*cylinderD, 2*cylinderD-fineExtend2[1]/2};
 
   auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
@@ -248,29 +279,13 @@ int main(int argc, char* argv[])
     fineGrid2.getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition2(fineGrid2.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition2);
-
-  prepareLattice(fineGrid, fineBulkDynamics, fineSBoundaryCondition);
-
-  const Vector<T,2> fineExtend3 {1.5*cylinderD, 1.5*cylinderD};
-  const Vector<T,2> fineOrigin3 {2*cylinderD-fineExtend3[0]/2, 2*cylinderD-fineExtend3[1]/2};
-
-  auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
-  prepareGeometry(fineGrid3);
-
-  disableRefinedArea(fineGrid2, fineGrid3);
-
-  BGKdynamics<T, DESCRIPTOR> fineBulkDynamics3(
-    fineGrid3.getConverter().getLatticeRelaxationFrequency(),
-    instances::getBulkMomenta<T, DESCRIPTOR>());
-
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition3(fineGrid3.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition3);
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineBoundaryCondition2(fineGrid2.getSuperLattice());
+  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineBoundaryCondition2);
 
-  prepareLattice(fineGrid2, fineBulkDynamics2, fineSBoundaryCondition2);
+  prepareLattice(fineGrid, fineBulkDynamics, fineBoundaryCondition);
+  prepareLattice(fineOutflowGrid, fineOutflowBulkDynamics, fineOutflowBoundaryCondition);
 
-  prepareLattice(fineGrid3, fineBulkDynamics3, fineSBoundaryCondition3);
+  prepareLattice(fineGrid2, fineBulkDynamics2, fineBoundaryCondition2);
 
   clout << "Starting simulation..." << endl;
   Timer<T> timer(
@@ -286,8 +301,8 @@ int main(int argc, char* argv[])
 
     getResults("level0_", coarseGrid, iT);
     getResults("level1_", fineGrid,   iT);
+    getResults("level1_outflow_", fineOutflowGrid, iT);
     getResults("level2_", fineGrid2,  iT);
-    getResults("level3_", fineGrid3,  iT);
 
     if (iT%statIter == 0) {
       timer.update(iT);