From 3edd576b152ec5df907597fb1a24d3e570d31f72 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Fri, 22 Mar 2019 22:25:48 +0100
Subject: Slightly shift cylinder relative to finest grid

This removes the "Bogus distance" warnings that are otherwise generated
by `addZeroVelocityBoundary`.
---
 apps/adrian/cylinder2d/common/model.h | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

(limited to 'apps/adrian/cylinder2d/common')

diff --git a/apps/adrian/cylinder2d/common/model.h b/apps/adrian/cylinder2d/common/model.h
index 124a1df..02e15e8 100644
--- a/apps/adrian/cylinder2d/common/model.h
+++ b/apps/adrian/cylinder2d/common/model.h
@@ -30,7 +30,7 @@ const T deltaR = cylinderD / N; // coarse lattice spacing
 const Vector<T,2> modelOrigin {0.0, 0.0};
 const Vector<T,2> modelExtend {22*cylinderD + deltaR, 4.1*cylinderD + deltaR};
 
-const Vector<T,2> cylinderCenter {2*cylinderD, 2*cylinderD + deltaR/2};
+const Vector<T,2> cylinderCenter {2*cylinderD, 2*cylinderD};
 
 void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid, Vector<T,2> origin, Vector<T,2> extend)
 {
@@ -78,7 +78,7 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid, Vector<T,2> origin, Vector<T,2>
 
   // Set material number for vertically centered cylinder
   {
-    const Vector<T,2> cylinderOrigin = origin + Vector<T,2> {cylinderCenter[0], cylinderCenter[1]};
+    const Vector<T,2> cylinderOrigin = origin + Vector<T,2> {cylinderCenter[0], cylinderCenter[1]+physSpacing/2};
     IndicatorCircle2D<T> obstacle(cylinderOrigin, cylinderD/2);
     sGeometry.rename(1,5,obstacle);
   }
@@ -99,6 +99,8 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid)
   auto& sGeometry = grid.getSuperGeometry();
   auto& sLattice  = grid.getSuperLattice();
 
+  const T physSpacing = converter.getPhysDeltaX();
+
   Dynamics<T,DESCRIPTOR>& bulkDynamics = grid.addDynamics(
       std::unique_ptr<Dynamics<T,DESCRIPTOR>>(
         new BGKdynamics<T,DESCRIPTOR>(
@@ -108,31 +110,26 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid)
   sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition = grid.getOnLatticeBoundaryCondition();
   createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);
 
-  const T omega = converter.getLatticeRelaxationFrequency();
+  sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& sOffBoundaryCondition = grid.getOffLatticeBoundaryCondition();
+  createBouzidiBoundaryCondition2D<T,DESCRIPTOR>(sOffBoundaryCondition);
 
   sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>());
   sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
   sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
   sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
   sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
+  sLattice.defineDynamics(sGeometry, 5, &instances::getNoDynamics<T,DESCRIPTOR>());
 
-  // only set zero velocity boundary if cylinder covered by the grid
-  if (sGeometry.getStatistics().getNvoxel(5) > 0) {
-    sLattice.defineDynamics(sGeometry, 5, &instances::getNoDynamics<T,DESCRIPTOR>());
-
-    sOffLatticeBoundaryCondition2D<T,DESCRIPTOR> sOffBoundaryCondition = grid.getOffLatticeBoundaryCondition();
-    createBouzidiBoundaryCondition2D<T,DESCRIPTOR>(sOffBoundaryCondition);
-
-    const Vector<T,2> cylinderOrigin {cylinderCenter[0], cylinderCenter[1]};
-    IndicatorCircle2D<T> obstacle(cylinderOrigin, cylinderD/2);
-
-    sOffBoundaryCondition.addZeroVelocityBoundary(sGeometry, 5, obstacle);
-  }
+  const T omega = converter.getLatticeRelaxationFrequency();
 
   sBoundaryCondition.addVelocityBoundary(sGeometry, 2, omega);
   sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega);
   sBoundaryCondition.addPressureBoundary(sGeometry, 4, omega);
 
+  const Vector<T,2> cylinderOrigin {cylinderCenter[0], cylinderCenter[1]+physSpacing/2};
+  IndicatorCircle2D<T> obstacle(cylinderOrigin, cylinderD/2);
+  sOffBoundaryCondition.addZeroVelocityBoundary(sGeometry, 5, obstacle);
+
   AnalyticalConst2D<T,T> rho0(1.0);
   AnalyticalConst2D<T,T> u0(0.0, 0.0);
 
@@ -143,6 +140,7 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid)
   sLattice.initialize();
 
   clout << "Prepare lattice ... OK" << std::endl;
+  sGeometry.print();
 }
 
 void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
@@ -201,14 +199,16 @@ void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid, int iT, bool print=true)
   auto& sGeometry = grid.getSuperGeometry();
   auto& converter = grid.getConverter();
 
+  const T physSpacing = converter.getPhysDeltaX();
+
   SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(sLattice, converter);
   AnalyticalFfromSuperF2D<T> intpolatePressure(pressure, true);
   SuperLatticePhysDrag2D<T,DESCRIPTOR> dragF(sLattice, sGeometry, 5, converter);
 
   const T radiusCylinder  = cylinderD/2;
 
-  const T point1[2] { cylinderCenter[0] - radiusCylinder, cylinderCenter[1] };
-  const T point2[2] { cylinderCenter[0] + radiusCylinder, cylinderCenter[1] };
+  const T point1[2] { cylinderCenter[0] - radiusCylinder, cylinderCenter[1]+physSpacing/2 };
+  const T point2[2] { cylinderCenter[0] + radiusCylinder, cylinderCenter[1]+physSpacing/2 };
 
   T pressureInFrontOfCylinder, pressureBehindCylinder;
   intpolatePressure(&pressureInFrontOfCylinder, point1);
-- 
cgit v1.2.3