Tidy up refined cylinder2d

1 files changed, 65 insertions, 54 deletions

diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp
index 2e18ebd..ffab6a0 100644
--- a/apps/adrian/cylinder2d/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/cylinder2d.cpp
@@ -31,21 +31,20 @@
 #include <vector>
 
 using namespace olb;
-using namespace olb::descriptors;
 
 typedef double T;
 
-#define DESCRIPTOR D2Q9Descriptor
+#define DESCRIPTOR descriptors::D2Q9Descriptor
 
 /// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
 const T cylinderD = 1.0;
 
-const int N = 10;              // resolution of the model
+const int N = 20;              // 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.00333;        // lattice velocity (also used as inflow velocity)
-const T maxPhysT = 120.;       // max. simulation time in s, SI unit
+const T uLat = 0.1;            // lattice velocity
+const T maxPhysT = 60.;        // max. simulation time in s, SI unit
 
 
 void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
@@ -88,7 +87,7 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
   {
     const Vector<T,2> origin {2*cylinderD, 2*cylinderD};
     IndicatorCircle2D<T> obstacle(origin, cylinderD/2);
-    sGeometry.rename(1,2,obstacle);
+    sGeometry.rename(1,5,obstacle);
   }
 
   sGeometry.clean();
@@ -99,6 +98,16 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
   clout << "Prepare Geometry ... OK" << std::endl;
 }
 
+void disableRefinedArea(Grid2D<T,DESCRIPTOR>& coarseGrid,
+                        RefiningGrid2D<T,DESCRIPTOR>& fineGrid)
+{
+  auto& sGeometry = coarseGrid.getSuperGeometry();
+  auto refinedOverlap = fineGrid.getRefinedOverlap();
+  sGeometry.rename(1,0,*refinedOverlap);
+  sGeometry.rename(2,0,*refinedOverlap);
+  sGeometry.rename(5,0,*refinedOverlap);
+}
+
 void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
                     Dynamics<T, DESCRIPTOR>& bulkDynamics,
                     sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition)
@@ -114,9 +123,10 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
 
   sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>());
   sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
-  sLattice.defineDynamics(sGeometry, 2, &instances::getNoDynamics<T,DESCRIPTOR>());
+  sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
   sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
   sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
+  sLattice.defineDynamics(sGeometry, 5, &instances::getBounceBack<T,DESCRIPTOR>()); // cylinder
 
   sBoundaryCondition.addVelocityBoundary(sGeometry, 2, omega);
   sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega);
@@ -125,7 +135,7 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
   AnalyticalConst2D<T,T> rho0(1.0);
   AnalyticalConst2D<T,T> u0(0.0, 0.0);
 
-  auto materials = sGeometry.getMaterialIndicator({1, 2, 3, 4});
+  auto materials = sGeometry.getMaterialIndicator({1, 2, 3, 4, 5});
   sLattice.defineRhoU(materials, rho0, u0);
   sLattice.iniEquilibrium(materials, rho0, u0);
 
@@ -140,11 +150,11 @@ void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
   auto& sGeometry = grid.getSuperGeometry();
   auto& sLattice  = grid.getSuperLattice();
 
-  const int iTmaxStart = converter.getLatticeTime(0.05 * maxPhysT);
-  const int iTupdate = 20;
+  const int iTmaxStart = converter.getLatticeTime(0.2*maxPhysT);
+  const int iTupdate = 10;
 
   if ( iT % iTupdate == 0 && iT <= iTmaxStart ) {
-    PolynomialStartScale<T,T> StartScale(iTmaxStart, 1.0);
+    PolynomialStartScale<T,T> StartScale(iTmaxStart, 1);
 
     T iTvec[1] { T(iT) };
     T frac[1] { };
@@ -162,7 +172,7 @@ void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
 
 void getResults(const std::string& prefix,
                 Grid2D<T,DESCRIPTOR>& grid,
-                int iT, Timer<T>& timer)
+                int iT)
 {
   OstreamManager clout(std::cout,"getResults");
 
@@ -180,8 +190,6 @@ void getResults(const std::string& prefix,
   vtmWriter.addFunctor(pressure);
   vtmWriter.addFunctor(knudsen);
 
-  const int statIter = converter.getLatticeTime(maxPhysT/10.);
-
   if (iT==0) {
     vtmWriter.createMasterFile();
   }
@@ -189,13 +197,6 @@ void getResults(const std::string& prefix,
   if (iT%100==0) {
     vtmWriter.write(iT);
   }
-
-  if (iT%statIter==0 ) {
-    timer.update(iT);
-    timer.printStep();
-
-    sLattice.getStatistics().print(iT,converter.getPhysTime(iT));
-  }
 }
 
 int main(int argc, char* argv[])
@@ -211,16 +212,6 @@ int main(int argc, char* argv[])
   Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, LatticeVelocity<T>(uLat), Re);
   prepareGeometry(coarseGrid);
 
-  const Vector<T,2> fineExtend {10*cylinderD, 2.75*cylinderD};
-  const Vector<T,2> fineOrigin {0.75*cylinderD, (ly-fineExtend[1])/2};
-
-  auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
-  prepareGeometry(fineGrid);
-
-  auto refinedOverlap = fineGrid.getRefinedOverlap();
-  coarseGrid.getSuperGeometry().rename(1,0,*refinedOverlap);
-  coarseGrid.getSuperGeometry().rename(2,0,*refinedOverlap);
-
   BGKdynamics<T, DESCRIPTOR> coarseBulkDynamics(
     coarseGrid.getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
@@ -228,7 +219,13 @@ int main(int argc, char* argv[])
   sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
   createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);
 
-  prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition);
+  const Vector<T,2> fineExtend {16*cylinderD, 3.75*cylinderD};
+  const Vector<T,2> fineOrigin {0.75*cylinderD, (ly-fineExtend[1])/2};
+
+  auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
+  prepareGeometry(fineGrid);
+
+  disableRefinedArea(coarseGrid, fineGrid);
 
   BGKdynamics<T, DESCRIPTOR> fineBulkDynamics(
     fineGrid.getConverter().getLatticeRelaxationFrequency(),
@@ -237,17 +234,15 @@ int main(int argc, char* argv[])
   sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
   createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
 
-  const Vector<T,2> fineExtend2 {1.75*cylinderD, 1.75*cylinderD};
-  const Vector<T,2> fineOrigin2 {2*cylinderD-fineExtend2[0]/2, 2*cylinderD-fineExtend2[1]/2};
+  prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition);
+
+  const Vector<T,2> fineExtend2 {11*cylinderD, 3*cylinderD};
+  const Vector<T,2> fineOrigin2 {1*cylinderD, 2*cylinderD-fineExtend2[1]/2};
 
   auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
   prepareGeometry(fineGrid2);
 
-  auto refinedOverlap2 = fineGrid2.getRefinedOverlap();
-  fineGrid.getSuperGeometry().rename(1,0,*refinedOverlap2);
-  fineGrid.getSuperGeometry().rename(2,0,*refinedOverlap2);
-
-  prepareLattice(fineGrid, fineBulkDynamics, fineSBoundaryCondition);
+  disableRefinedArea(fineGrid, fineGrid2);
 
   BGKdynamics<T, DESCRIPTOR> fineBulkDynamics2(
     fineGrid2.getConverter().getLatticeRelaxationFrequency(),
@@ -256,34 +251,50 @@ int main(int argc, char* argv[])
   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);
+
   prepareLattice(fineGrid2, fineBulkDynamics2, fineSBoundaryCondition2);
 
+  prepareLattice(fineGrid3, fineBulkDynamics3, fineSBoundaryCondition3);
+
   clout << "Starting simulation..." << endl;
   Timer<T> timer(
     coarseGrid.getConverter().getLatticeTime(maxPhysT),
     coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
   timer.start();
+  const int statIter = coarseGrid.getConverter().getLatticeTime(0.5);
 
   for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
     setBoundaryValues(coarseGrid, iT);
 
     coarseGrid.collideAndStream();
 
-    getResults(
-      "coarse_",
-      coarseGrid,
-      iT,
-      timer);
-    getResults(
-      "fine_",
-      fineGrid,
-      iT,
-      timer);
-    getResults(
-      "fine2_",
-      fineGrid2,
-      iT,
-      timer);
+    getResults("level0_", coarseGrid, iT);
+    getResults("level1_", fineGrid,   iT);
+    getResults("level2_", fineGrid2,  iT);
+    getResults("level3_", fineGrid3,  iT);
+
+    if (iT%statIter == 0) {
+      timer.update(iT);
+      timer.printStep();
+
+      coarseGrid.getSuperLattice().getStatistics().print(iT, coarseGrid.getConverter().getPhysTime(iT));
+    }
   }
 
   timer.stop();