summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAdrian Kummerlaender2019-02-03 19:51:42 +0100
committerAdrian Kummerlaender2019-06-24 15:17:28 +0200
commit953db99dca6a34c0e9e6f71561f8285b08822bcc (patch)
tree2d1cde3c00c23e1195eb92ddd81804d8d2731d30
parent12cde59008e15abb935a2177a2a4ed3038841c96 (diff)
downloadgrid_refinement_openlb-953db99dca6a34c0e9e6f71561f8285b08822bcc.tar
grid_refinement_openlb-953db99dca6a34c0e9e6f71561f8285b08822bcc.tar.gz
grid_refinement_openlb-953db99dca6a34c0e9e6f71561f8285b08822bcc.tar.bz2
grid_refinement_openlb-953db99dca6a34c0e9e6f71561f8285b08822bcc.tar.xz
grid_refinement_openlb-953db99dca6a34c0e9e6f71561f8285b08822bcc.zip
Revamp parametrization of refined cylinder2d geometry
Finally seems to stop mixing up material numbers for every other resolution I try it out with. Cylinder diameter is now actually set to 0.1m as called for by [SchaeferTurek96].
-rw-r--r--apps/adrian/cylinder2d/cylinder2d.cpp178
-rw-r--r--apps/adrian/poiseuille2d/poiseuille2d.cpp2
2 files changed, 94 insertions, 86 deletions
diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp
index c853939..066bb6f 100644
--- a/apps/adrian/cylinder2d/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/cylinder2d.cpp
@@ -37,17 +37,25 @@ typedef double T;
#define DESCRIPTOR descriptors::D2Q9Descriptor
/// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
-const T cylinderD = 1.0;
-
-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.05; // lattice velocity
-const T maxPhysT = 60.; // max. simulation time in s, SI unit
-
-
-void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
+const T cylinderD = 0.1;
+const int N = 5; // 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 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. density
+
+void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid, Vector<T,2> origin, Vector<T,2> extend)
{
OstreamManager clout(std::cout,"prepareGeometry");
clout << "Prepare Geometry ..." << std::endl;
@@ -59,34 +67,42 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
const T physSpacing = converter.getPhysDeltaX();
- // Set material number for bounce back boundaries
+ // Set material number for channel walls
{
- const Vector<T,2> wallExtend {lx+physSpacing, physSpacing};
- const Vector<T,2> wallOrigin {-physSpacing/2, -physSpacing/2};
+ const Vector<T,2> wallExtend { extend[0]+physSpacing/2, physSpacing/2 };
+ const Vector<T,2> wallOrigin = origin - physSpacing/4;
IndicatorCuboid2D<T> lowerWall(wallExtend, wallOrigin);
sGeometry.rename(1,2,lowerWall);
+ }
+ {
+ const Vector<T,2> wallExtend { extend[0]+physSpacing/2, physSpacing/2 };
+ const Vector<T,2> wallOrigin { origin[0]-physSpacing/4, extend[1]-physSpacing/4 };
- IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin + Vector<T,2> {0,ly-physSpacing/2});
+ IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin);
sGeometry.rename(1,2,upperWall);
}
// Set material number for inflow and outflow
{
- const Vector<T,2> extend { physSpacing/2, ly-physSpacing/2};
- const Vector<T,2> origin {-physSpacing/4, -physSpacing/4};
+ const Vector<T,2> inflowExtend { physSpacing/2, extend[1]+physSpacing/4 };
+ const Vector<T,2> inflowOrigin = origin - physSpacing/4;
- IndicatorCuboid2D<T> inflow(extend, origin);
+ IndicatorCuboid2D<T> inflow(inflowExtend, inflowOrigin);
sGeometry.rename(1,3,inflow);
+ }
+ {
+ const Vector<T,2> outflowExtend { physSpacing/2, extend[1]+physSpacing/4 };
+ const Vector<T,2> outflowOrigin { extend[0]-physSpacing/4, origin[0]-physSpacing/4 };
- IndicatorCuboid2D<T> outflow(extend, origin + Vector<T,2> {lx,0});
+ IndicatorCuboid2D<T> outflow(outflowExtend, outflowOrigin);
sGeometry.rename(1,4,outflow);
}
// Set material number for vertically centered cylinder
{
- const Vector<T,2> origin {2*cylinderD, 2*cylinderD};
- IndicatorCircle2D<T> obstacle(origin, cylinderD/2);
+ const Vector<T,2> cylinderOrigin = origin + Vector<T,2> {cylinderX, cylinderY};
+ IndicatorCircle2D<T> obstacle(cylinderOrigin, cylinderD/2);
sGeometry.rename(1,5,obstacle);
}
@@ -121,25 +137,25 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid)
instances::getBulkMomenta<T,DESCRIPTOR>())));
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition = grid.getOnLatticeBoundaryCondition();
- createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);
+ createInterpBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);
const T omega = converter.getLatticeRelaxationFrequency();
sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>());
sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
- sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
+ sLattice.defineDynamics(sGeometry, 2, &instances::getBounceBack<T,DESCRIPTOR>()); // walls
sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
- sLattice.defineDynamics(sGeometry, 5, &instances::getBounceBack<T,DESCRIPTOR>()); // cylinder
+ sLattice.defineDynamics(sGeometry, 5, &bulkDynamics); // 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, 2, 3, 4, 5});
+ auto materials = sGeometry.getMaterialIndicator({1, 3, 4, 5});
sLattice.defineRhoU(materials, rho0, u0);
sLattice.iniEquilibrium(materials, rho0, u0);
@@ -155,8 +171,8 @@ void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
auto& sGeometry = grid.getSuperGeometry();
auto& sLattice = grid.getSuperLattice();
- const int iTmaxStart = converter.getLatticeTime(0.2*maxPhysT);
- const int iTupdate = 10;
+ const int iTmaxStart = converter.getLatticeTime(0.4*16);
+ const int iTupdate = 5;
if ( iT % iTupdate == 0 && iT <= iTmaxStart ) {
PolynomialStartScale<T,T> StartScale(iTmaxStart, 1);
@@ -165,11 +181,8 @@ void setBoundaryValues(Grid2D<T,DESCRIPTOR>& grid, int iT)
T frac[1] { };
StartScale(frac, iTvec);
- const T maxVelocity = converter.getCharLatticeVelocity() * frac[0];
- const T radius = ly/2;
- std::vector<T> axisPoint{0, ly/2};
- std::vector<T> axisDirection{1, 0};
- Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);
+ const T maxVelocity = converter.getCharLatticeVelocity() * 3./2. * frac[0];
+ Poiseuille2D<T> u(sGeometry, 3, maxVelocity, deltaR/2);
sLattice.defineU(sGeometry, 3, u);
}
@@ -199,13 +212,13 @@ void getResults(const std::string& prefix,
vtmWriter.createMasterFile();
}
- if (iT%100==0) {
- vtmWriter.write(iT);
- }
+ vtmWriter.write(iT);
}
void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid)
{
+ static T maxDrag = 0.0;
+
OstreamManager clout(std::cout,"measurement");
auto& sLattice = grid.getSuperLattice();
@@ -214,14 +227,12 @@ void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid)
SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(sLattice, converter);
AnalyticalFfromSuperF2D<T> intpolatePressure(pressure, true);
- SuperLatticePhysDrag2D<T,DESCRIPTOR> drag(sLattice, sGeometry, 5, converter);
+ SuperLatticePhysDrag2D<T,DESCRIPTOR> dragF(sLattice, sGeometry, 5, converter);
- const T centerCylinderX = 2*cylinderD;
- const T centerCylinderY = 2*cylinderD;
const T radiusCylinder = cylinderD/2;
- const T point1[2] { centerCylinderX - radiusCylinder, centerCylinderY };
- const T point2[2] { centerCylinderX + radiusCylinder, centerCylinderY };
+ const T point1[2] { cylinderX - radiusCylinder, cylinderY };
+ const T point2[2] { cylinderX + radiusCylinder, cylinderY };
T pressureInFrontOfCylinder, pressureBehindCylinder;
intpolatePressure(&pressureInFrontOfCylinder, point1);
@@ -234,9 +245,12 @@ void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid)
clout << "; pressureDrop=" << pressureDrop;
const int input[3] {};
- T _drag[drag.getTargetDim()] {};
- drag(_drag, input);
- clout << "; drag=" << _drag[0] << "; lift=" << _drag[1] << endl;
+ T drag[dragF.getTargetDim()] {};
+ dragF(drag, input);
+ if (drag[0] > maxDrag) {
+ maxDrag = drag[0];
+ };
+ clout << "; drag=" << drag[0] << "; maxDrag: " << maxDrag << "; lift=" << drag[1] << endl;
}
int main(int argc, char* argv[])
@@ -249,83 +263,77 @@ int main(int argc, char* argv[])
const Vector<T,2> coarseExtend {lx, ly};
IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
- Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, LatticeVelocity<T>(uLat), Re);
- prepareGeometry(coarseGrid);
+ Grid2D<T,DESCRIPTOR> coarseGrid(
+ coarseCuboid,
+ RelaxationTime<T>(tau),
+ N,
+ PhysCharacteristics);
+ const Vector<T,2> domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0);
+ const Vector<T,2> domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0);
+
+ prepareGeometry(coarseGrid, domainOrigin, domainExtend);
const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
- const Vector<T,2> fineExtend {6.5*cylinderD, ly-2*coarseDeltaX};
- const Vector<T,2> fineOrigin {0.75*cylinderD, coarseDeltaX};
+ const Vector<T,2> fineExtend {10*cylinderD, ly-2*coarseDeltaX};
+ const Vector<T,2> fineOrigin {1*coarseDeltaX, coarseDeltaX};
auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
- prepareGeometry(fineGrid);
+ prepareGeometry(fineGrid, domainOrigin, domainExtend);
disableRefinedArea(coarseGrid, fineGrid);
- const Vector<T,2> fineOutflowExtend {5*coarseDeltaX, ly};
- const Vector<T,2> fineOutflowOrigin {lx-5*coarseDeltaX, 0};
-
- auto& fineOutflowGrid = coarseGrid.refine(fineOutflowOrigin, fineOutflowExtend, 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().reset(refined);
- }
-
- const Vector<T,2> fineExtend2 {2.3*cylinderD, 1.7*cylinderD};
- const Vector<T,2> fineOrigin2 {1.0*cylinderD, 2*cylinderD-fineExtend2[1]/2};
+ const Vector<T,2> fineExtend2 {6*cylinderD, fineGrid.getExtend()[1]-2*coarseDeltaX};
+ const Vector<T,2> fineOrigin2 {2*coarseDeltaX, (ly-fineExtend2[1])/2};
auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
- prepareGeometry(fineGrid2);
+ prepareGeometry(fineGrid2, domainOrigin, domainExtend);
disableRefinedArea(fineGrid, fineGrid2);
- const Vector<T,2> fineExtend3 {1.4*cylinderD, 1.4*cylinderD};
- const Vector<T,2> fineOrigin3 {2*cylinderD-fineExtend3[0]/2, 2*cylinderD-fineExtend3[1]/2};
+ const Vector<T,2> fineExtend3 {4*cylinderD, 2*cylinderD};
+ const Vector<T,2> fineOrigin3 {1.2*cylinderD, (ly-fineExtend3[1])/2};
auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
- prepareGeometry(fineGrid3);
+ prepareGeometry(fineGrid3, domainOrigin, domainExtend);
disableRefinedArea(fineGrid2, fineGrid3);
+ const Vector<T,2> fineExtend4 {1.25*cylinderD, 1.25*cylinderD};
+ const Vector<T,2> fineOrigin4 {cylinderX - fineExtend4[0]/2, cylinderY - fineExtend4[1]/2};
+
+ auto& fineGrid4 = fineGrid3.refine(fineOrigin4, fineExtend4);
+ prepareGeometry(fineGrid4, domainOrigin, domainExtend);
+ disableRefinedArea(fineGrid3, fineGrid4);
+
prepareLattice(coarseGrid);
prepareLattice(fineGrid);
- prepareLattice(fineOutflowGrid);
prepareLattice(fineGrid2);
prepareLattice(fineGrid3);
+ prepareLattice(fineGrid4);
clout << "Total number of active cells: " << coarseGrid.getActiveVoxelN() << endl;
clout << "Starting simulation..." << endl;
+
+ const int statIter = coarseGrid.getConverter().getLatticeTime(0.1);
Timer<T> timer(
coarseGrid.getConverter().getLatticeTime(maxPhysT),
coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
timer.start();
- const int statIter = coarseGrid.getConverter().getLatticeTime(1);
- for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
+ for (int iT = 0; iT <= coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
setBoundaryValues(coarseGrid, iT);
coarseGrid.collideAndStream();
- 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);
timer.printStep();
- coarseGrid.getSuperLattice().getStatistics().print(iT, coarseGrid.getConverter().getPhysTime(iT));
+ getResults("level0_", coarseGrid, iT);
+ getResults("level1_", fineGrid, iT);
+ getResults("level2_", fineGrid2, iT);
+ getResults("level3_", fineGrid3, iT);
+ getResults("level4_", fineGrid4, iT);
- takeMeasurements(fineGrid3);
+ takeMeasurements(fineGrid4);
}
}
diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 5a47f96..d4d8e0f 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -222,7 +222,7 @@ int main(int argc, char* argv[])
const Vector<T,2> coarseExtend {lx/2, ly};
IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
- Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, LatticeVelocity<T>(uMax), Re);
+ Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, LatticeVelocity<T>(uMax), N, Re);
prepareGeometry(coarseGrid);
const T coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();