Replicate SchaeferTurek 2D-2 setup in refined cylinder2d

Weirdly it seems easier to get reasonably close values for the unsteady
flow. Results for the steady flow seem to get worse when increasing
grid resolution even in uniform grids.

1 files changed, 15 insertions, 18 deletions

diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp
index 066bb6f..5b67408 100644
--- a/apps/adrian/cylinder2d/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/cylinder2d.cpp
@@ -38,21 +38,21 @@ typedef double T;
 
 /// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
 const T   cylinderD = 0.1;
-const int N         = 5;                      // resolution of the cylinder
+const int N         = 10;                      // resolution of the cylinder
 const T   deltaR    = cylinderD / N;          // coarse lattice spacing
 const T   lx        = 22*cylinderD + deltaR;  // length of the channel
 const T   ly        = 4.1*cylinderD + deltaR; // height of the channel
 const T   cylinderX = 2*cylinderD;
 const T   cylinderY = 2*cylinderD + deltaR/2;
 
-const T   Re        = 20.;                    // Reynolds number
+const T   Re        = 100.;                   // Reynolds number
 const T   tau       = 0.51;                   // relaxation time
 const T   maxPhysT  = 16.;                    // max. simulation time in s, SI unit
 
 const Characteristics<T> PhysCharacteristics(
   0.1,         // char. phys. length
-  0.2,         // char. phys. velocity
-  0.1*0.2/Re,  // phsy. kinematic viscosity
+  1.0,         // char. phys. velocity
+  0.1/Re,      // phsy. kinematic viscosity
   1.0);        // char. phys. density
 
 void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid, Vector<T,2> origin, Vector<T,2> extend)
@@ -143,19 +143,19 @@ 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::getBounceBack<T,DESCRIPTOR>()); // walls
+  sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
   sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
   sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
-  sLattice.defineDynamics(sGeometry, 5, &bulkDynamics); // cylinder
+  sLattice.defineDynamics(sGeometry, 5, &instances::getBounceBack<T,DESCRIPTOR>()); // cylinder
 
+  sBoundaryCondition.addVelocityBoundary(sGeometry, 2, omega);
   sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega);
   sBoundaryCondition.addPressureBoundary(sGeometry, 4, omega);
-  sBoundaryCondition.addVelocityBoundary(sGeometry, 5, omega);
 
   AnalyticalConst2D<T,T> rho0(1.0);
   AnalyticalConst2D<T,T> u0(0.0, 0.0);
 
-  auto materials = sGeometry.getMaterialIndicator({1, 3, 4, 5});
+  auto materials = sGeometry.getMaterialIndicator({1, 2, 3, 4});
   sLattice.defineRhoU(materials, rho0, u0);
   sLattice.iniEquilibrium(materials, rho0, u0);
 
@@ -238,11 +238,8 @@ void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid)
   intpolatePressure(&pressureInFrontOfCylinder, point1);
   intpolatePressure(&pressureBehindCylinder,    point2);
 
-  clout << "pressure1=" << pressureInFrontOfCylinder;
-  clout << "; pressure2=" << pressureBehindCylinder;
-
   T pressureDrop = pressureInFrontOfCylinder - pressureBehindCylinder;
-  clout << "; pressureDrop=" << pressureDrop;
+  clout << "pressureDrop=" << pressureDrop;
 
   const int input[3] {};
   T drag[dragF.getTargetDim()] {};
@@ -275,22 +272,22 @@ int main(int argc, char* argv[])
 
   const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
 
-  const Vector<T,2> fineExtend {10*cylinderD, ly-2*coarseDeltaX};
-  const Vector<T,2> fineOrigin {1*coarseDeltaX, coarseDeltaX};
+  const Vector<T,2> fineExtend {10*cylinderD, domainExtend[1]-4*coarseDeltaX};
+  const Vector<T,2> fineOrigin {0.5*cylinderD, (domainExtend[1]-fineExtend[1])/2};
 
   auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
   prepareGeometry(fineGrid, domainOrigin, domainExtend);
   disableRefinedArea(coarseGrid, fineGrid);
 
-  const Vector<T,2> fineExtend2 {6*cylinderD, fineGrid.getExtend()[1]-2*coarseDeltaX};
-  const Vector<T,2> fineOrigin2 {2*coarseDeltaX, (ly-fineExtend2[1])/2};
+  const Vector<T,2> fineExtend2 {6*cylinderD, fineGrid.getExtend()[1]-4*coarseDeltaX};
+  const Vector<T,2> fineOrigin2 {0.75*cylinderD, (domainExtend[1]-fineExtend2[1])/2};
 
   auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
   prepareGeometry(fineGrid2, domainOrigin, domainExtend);
   disableRefinedArea(fineGrid, fineGrid2);
 
   const Vector<T,2> fineExtend3 {4*cylinderD, 2*cylinderD};
-  const Vector<T,2> fineOrigin3 {1.2*cylinderD, (ly-fineExtend3[1])/2};
+  const Vector<T,2> fineOrigin3 {1*cylinderD, (domainExtend[1]-fineExtend3[1])/2};
 
   auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
   prepareGeometry(fineGrid3, domainOrigin, domainExtend);
@@ -312,7 +309,7 @@ int main(int argc, char* argv[])
   clout << "Total number of active cells: " << coarseGrid.getActiveVoxelN() << endl;
   clout << "Starting simulation..." << endl;
 
-  const int statIter = coarseGrid.getConverter().getLatticeTime(0.1);
+  const int statIter = coarseGrid.getConverter().getLatticeTime(0.05);
   Timer<T> timer(
     coarseGrid.getConverter().getLatticeTime(maxPhysT),
     coarseGrid.getSuperGeometry().getStatistics().getNvoxel());