From 953db99dca6a34c0e9e6f71561f8285b08822bcc Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sun, 3 Feb 2019 19:51:42 +0100 Subject: 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]. --- apps/adrian/cylinder2d/cylinder2d.cpp | 178 ++++++++++++++++-------------- apps/adrian/poiseuille2d/poiseuille2d.cpp | 2 +- 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& 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 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& grid, Vector origin, Vector extend) { OstreamManager clout(std::cout,"prepareGeometry"); clout << "Prepare Geometry ..." << std::endl; @@ -59,34 +67,42 @@ void prepareGeometry(Grid2D& grid) const T physSpacing = converter.getPhysDeltaX(); - // Set material number for bounce back boundaries + // Set material number for channel walls { - const Vector wallExtend {lx+physSpacing, physSpacing}; - const Vector wallOrigin {-physSpacing/2, -physSpacing/2}; + const Vector wallExtend { extend[0]+physSpacing/2, physSpacing/2 }; + const Vector wallOrigin = origin - physSpacing/4; IndicatorCuboid2D lowerWall(wallExtend, wallOrigin); sGeometry.rename(1,2,lowerWall); + } + { + const Vector wallExtend { extend[0]+physSpacing/2, physSpacing/2 }; + const Vector wallOrigin { origin[0]-physSpacing/4, extend[1]-physSpacing/4 }; - IndicatorCuboid2D upperWall(wallExtend, wallOrigin + Vector {0,ly-physSpacing/2}); + IndicatorCuboid2D upperWall(wallExtend, wallOrigin); sGeometry.rename(1,2,upperWall); } // Set material number for inflow and outflow { - const Vector extend { physSpacing/2, ly-physSpacing/2}; - const Vector origin {-physSpacing/4, -physSpacing/4}; + const Vector inflowExtend { physSpacing/2, extend[1]+physSpacing/4 }; + const Vector inflowOrigin = origin - physSpacing/4; - IndicatorCuboid2D inflow(extend, origin); + IndicatorCuboid2D inflow(inflowExtend, inflowOrigin); sGeometry.rename(1,3,inflow); + } + { + const Vector outflowExtend { physSpacing/2, extend[1]+physSpacing/4 }; + const Vector outflowOrigin { extend[0]-physSpacing/4, origin[0]-physSpacing/4 }; - IndicatorCuboid2D outflow(extend, origin + Vector {lx,0}); + IndicatorCuboid2D outflow(outflowExtend, outflowOrigin); sGeometry.rename(1,4,outflow); } // Set material number for vertically centered cylinder { - const Vector origin {2*cylinderD, 2*cylinderD}; - IndicatorCircle2D obstacle(origin, cylinderD/2); + const Vector cylinderOrigin = origin + Vector {cylinderX, cylinderY}; + IndicatorCircle2D obstacle(cylinderOrigin, cylinderD/2); sGeometry.rename(1,5,obstacle); } @@ -121,25 +137,25 @@ void prepareLattice(Grid2D& grid) instances::getBulkMomenta()))); sOnLatticeBoundaryCondition2D& sBoundaryCondition = grid.getOnLatticeBoundaryCondition(); - createLocalBoundaryCondition2D(sBoundaryCondition); + createInterpBoundaryCondition2D(sBoundaryCondition); const T omega = converter.getLatticeRelaxationFrequency(); sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics()); sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk - sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls + sLattice.defineDynamics(sGeometry, 2, &instances::getBounceBack()); // walls sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow - sLattice.defineDynamics(sGeometry, 5, &instances::getBounceBack()); // 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 rho0(1.0); AnalyticalConst2D 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& 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 StartScale(iTmaxStart, 1); @@ -165,11 +181,8 @@ void setBoundaryValues(Grid2D& grid, int iT) T frac[1] { }; StartScale(frac, iTvec); - const T maxVelocity = converter.getCharLatticeVelocity() * frac[0]; - const T radius = ly/2; - std::vector axisPoint{0, ly/2}; - std::vector axisDirection{1, 0}; - Poiseuille2D u(axisPoint, axisDirection, maxVelocity, radius); + const T maxVelocity = converter.getCharLatticeVelocity() * 3./2. * frac[0]; + Poiseuille2D 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& grid) { + static T maxDrag = 0.0; + OstreamManager clout(std::cout,"measurement"); auto& sLattice = grid.getSuperLattice(); @@ -214,14 +227,12 @@ void takeMeasurements(Grid2D& grid) SuperLatticePhysPressure2D pressure(sLattice, converter); AnalyticalFfromSuperF2D intpolatePressure(pressure, true); - SuperLatticePhysDrag2D drag(sLattice, sGeometry, 5, converter); + SuperLatticePhysDrag2D 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& 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 coarseExtend {lx, ly}; IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin); - Grid2D coarseGrid(coarseCuboid, N, LatticeVelocity(uLat), Re); - prepareGeometry(coarseGrid); + Grid2D coarseGrid( + coarseCuboid, + RelaxationTime(tau), + N, + PhysCharacteristics); + const Vector domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0); + const Vector domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0); + + prepareGeometry(coarseGrid, domainOrigin, domainExtend); const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX(); - const Vector fineExtend {6.5*cylinderD, ly-2*coarseDeltaX}; - const Vector fineOrigin {0.75*cylinderD, coarseDeltaX}; + const Vector fineExtend {10*cylinderD, ly-2*coarseDeltaX}; + const Vector fineOrigin {1*coarseDeltaX, coarseDeltaX}; auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend); - prepareGeometry(fineGrid); + prepareGeometry(fineGrid, domainOrigin, domainExtend); disableRefinedArea(coarseGrid, fineGrid); - const Vector fineOutflowExtend {5*coarseDeltaX, ly}; - const Vector fineOutflowOrigin {lx-5*coarseDeltaX, 0}; - - auto& fineOutflowGrid = coarseGrid.refine(fineOutflowOrigin, fineOutflowExtend, false); - prepareGeometry(fineOutflowGrid); - - { - const Vector origin = fineOutflowGrid.getOrigin(); - const Vector extend = fineOutflowGrid.getExtend(); - - const Vector extendY = {0,extend[1]}; - - coarseGrid.addFineCoupling(fineOutflowGrid, origin, extendY); - coarseGrid.addCoarseCoupling(fineOutflowGrid, origin + Vector {coarseDeltaX,0}, extendY); - - IndicatorCuboid2D refined(extend, origin + Vector {2*coarseDeltaX,0}); - coarseGrid.getSuperGeometry().reset(refined); - } - - const Vector fineExtend2 {2.3*cylinderD, 1.7*cylinderD}; - const Vector fineOrigin2 {1.0*cylinderD, 2*cylinderD-fineExtend2[1]/2}; + const Vector fineExtend2 {6*cylinderD, fineGrid.getExtend()[1]-2*coarseDeltaX}; + const Vector fineOrigin2 {2*coarseDeltaX, (ly-fineExtend2[1])/2}; auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2); - prepareGeometry(fineGrid2); + prepareGeometry(fineGrid2, domainOrigin, domainExtend); disableRefinedArea(fineGrid, fineGrid2); - const Vector fineExtend3 {1.4*cylinderD, 1.4*cylinderD}; - const Vector fineOrigin3 {2*cylinderD-fineExtend3[0]/2, 2*cylinderD-fineExtend3[1]/2}; + const Vector fineExtend3 {4*cylinderD, 2*cylinderD}; + const Vector 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 fineExtend4 {1.25*cylinderD, 1.25*cylinderD}; + const Vector 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 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 coarseExtend {lx/2, ly}; IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin); - Grid2D coarseGrid(coarseCuboid, N, LatticeVelocity(uMax), Re); + Grid2D coarseGrid(coarseCuboid, LatticeVelocity(uMax), N, Re); prepareGeometry(coarseGrid); const T coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX(); -- cgit v1.2.3