Setup geometry relative to cylinder diameter

1 files changed, 25 insertions, 23 deletions

diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp
index 33028c0..3baa4ef 100644
--- a/apps/adrian/cylinder2d/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/cylinder2d.cpp
@@ -37,12 +37,16 @@ typedef double T;
 
 #define DESCRIPTOR D2Q9Descriptor
 
-const T lx  = 8.0;            // length of the channel
-const T ly  = 2.0;            // height of the channel
-const int N = 100;            // resolution of the model
-const T Re  = 1000.;          // Reynolds number
-const T baseTau = 0.52;       // Relaxation time of coarsest grid
-const T maxPhysT = 120.;      // max. simulation time in s, SI unit
+/// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
+const T cylinderD = 1.0;
+
+const int N = 10;              // 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
+
 
 void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
 {
@@ -58,8 +62,8 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
 
   // Set material number for bounce back boundaries
   {
-    const Vector<T,2> wallExtend {lx+physSpacing,  physSpacing/2};
-    const Vector<T,2> wallOrigin {-physSpacing/4, -physSpacing/4};
+    const Vector<T,2> wallExtend {lx+physSpacing,  physSpacing};
+    const Vector<T,2> wallOrigin {-physSpacing/2, -physSpacing/2};
 
     IndicatorCuboid2D<T> lowerWall(wallExtend, wallOrigin);
     sGeometry.rename(1,2,lowerWall);
@@ -82,8 +86,8 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
 
   // Set material number for vertically centered cylinder
   {
-    const Vector<T,2> origin {1.25, ly/2};
-    IndicatorCircle2D<T> obstacle(origin, 0.2);
+    const Vector<T,2> origin {2*cylinderD, 2*cylinderD};
+    IndicatorCircle2D<T> obstacle(origin, cylinderD/2);
     sGeometry.rename(1,2,obstacle);
   }
 
@@ -108,19 +112,20 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
 
   const T omega = converter.getLatticeRelaxationFrequency();
 
-  sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
+  sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>());
   sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
-  sLattice.defineDynamics(sGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
+  sLattice.defineDynamics(sGeometry, 2, &instances::getNoDynamics<T,DESCRIPTOR>());
   sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
   sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
 
+  sBoundaryCondition.addVelocityBoundary(sGeometry, 2, omega);
   sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega);
   sBoundaryCondition.addPressureBoundary(sGeometry, 4, omega);
 
   AnalyticalConst2D<T,T> rho0(1.0);
   AnalyticalConst2D<T,T> u0(0.0, 0.0);
 
-  auto materials = sGeometry.getMaterialIndicator({1, 3, 4});
+  auto materials = sGeometry.getMaterialIndicator({1, 2, 3, 4});
   sLattice.defineRhoU(materials, rho0, u0);
   sLattice.iniEquilibrium(materials, rho0, u0);
 
@@ -135,13 +140,10 @@ void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
   auto& sGeometry = grid.getSuperGeometry();
   auto& sLattice  = grid.getSuperLattice();
 
-  const int iTmaxStart = converter.getLatticeTime(0.1 * maxPhysT);
-  const int iTupdate = 50;
+  const int iTmaxStart = converter.getLatticeTime(0.05 * maxPhysT);
+  const int iTupdate = 20;
 
   if ( iT % iTupdate == 0 && iT <= iTmaxStart ) {
-    OstreamManager clout(std::cout,"setBoundaryValues");
-    clout << "Update inflow." << std::endl;
-
     PolynomialStartScale<T,T> StartScale(iTmaxStart, 1.0);
 
     T iTvec[1] { T(iT) };
@@ -204,11 +206,11 @@ int main(int argc, char* argv[])
   const Vector<T,2> coarseExtend {lx, ly};
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
 
-  Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, RelaxationTime<T>(baseTau), Re);
+  Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, LatticeVelocity<T>(uLat), Re);
   prepareGeometry(coarseGrid);
 
-  const Vector<T,2> fineExtend {3.0, 1.5};
-  const Vector<T,2> fineOrigin {0.8, (ly-fineExtend[1])/2};
+  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);
@@ -233,8 +235,8 @@ int main(int argc, char* argv[])
   sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
   createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
 
-  const Vector<T,2> fineExtend2 {1.2, 0.8};
-  const Vector<T,2> fineOrigin2 {0.95, (ly-fineExtend2[1])/2};
+  const Vector<T,2> fineExtend2 {5*cylinderD, 1.5*cylinderD};
+  const Vector<T,2> fineOrigin2 {1.25*cylinderD, (ly-fineExtend2[1])/2};
 
   auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
   prepareGeometry(fineGrid2);