From d0f1e440eaa92116781f69c953db47f758bde9b4 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Fri, 4 Jan 2019 18:02:54 +0100 Subject: Change F2C restriction, some cleanup --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 287 ++++++++++++++++-------------- 1 file changed, 150 insertions(+), 137 deletions(-) (limited to 'apps') diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index 8bb982d..22c35cb 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -27,25 +27,19 @@ #include "olb2D.hh" #include -#include -#include -#include -#include using namespace olb; using namespace olb::descriptors; -using namespace olb::graphics; -using namespace std; typedef double T; #define DESCRIPTOR D2Q9Descriptor -const T lx = 2.; // length of the channel +const T lx = 4.; // length of the channel const T ly = 1.; // height of the channel -const int N = 16; // resolution of the model -const T Re = 10.; // Reynolds number -const T maxPhysT = 40.; // max. simulation time in s, SI unit +const int N = 30; // resolution of the model +const T Re = 100.; // Reynolds number +const T maxPhysT = 60.; // max. simulation time in s, SI unit const T physInterval = 0.25; // interval for the convergence check in s const T residuum = 1e-5; // residuum for the convergence check @@ -61,14 +55,11 @@ void prepareGeometry(UnitConverter const& converter, // Set material number for bounce back boundaries superGeometry.rename(2,1,0,1); - Vector extend; - Vector origin; T physSpacing = converter.getPhysDeltaX(); + Vector extend{ physSpacing / 2, ly }; + Vector origin{ -physSpacing / 4, 0 }; // Set material number for inflow - extend[1] = ly; - extend[0] = physSpacing / 2; - origin[0] -= physSpacing / 4; IndicatorCuboid2D inflow(extend, origin); superGeometry.rename(1,3,1,inflow); @@ -77,6 +68,9 @@ void prepareGeometry(UnitConverter const& converter, IndicatorCuboid2D outflow(extend, origin); superGeometry.rename(1,4,1,outflow); + //IndicatorCuboid2D obstacle(Vector{0.5,0.2}, Vector{1.75,0.4}); + //superGeometry.rename(1,2,obstacle); + // Removes all not needed boundary voxels outside the surface superGeometry.clean(); // Removes all not needed boundary voxels inside the surface @@ -106,30 +100,28 @@ void prepareCoarseLattice(UnitConverter const& converter, sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega); - T Lx = converter.getLatticeLength(lx); - T Ly = converter.getLatticeLength(ly); - - T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); + const T Lx = converter.getLatticeLength(lx); + const T Ly = converter.getLatticeLength(ly); + const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1); - AnalyticalConst2D iniRho(1); - AnalyticalConst2D iniU(std::vector {0.,0.}); - - sLattice.defineRhoU(superGeometry, 1, iniRho, iniU); - sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU); - sLattice.defineRhoU(superGeometry, 2, iniRho, iniU); - sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU); + const T maxVelocity = converter.getCharLatticeVelocity(); + const T radius = ly/2; + std::vector axisPoint{0, ly/2}; + std::vector axisDirection{1, 0}; + Poiseuille2D u(axisPoint, axisDirection, maxVelocity, radius); - T maxVelocity = converter.getCharLatticeVelocity(); - T distance2Wall = converter.getConversionFactorLength(); - Poiseuille2D u(superGeometry, 3, maxVelocity, distance2Wall); + sLattice.defineRhoU(superGeometry, 1, rho, u); + sLattice.iniEquilibrium(superGeometry, 1, rho, u); + sLattice.defineRhoU(superGeometry, 2, rho, u); + sLattice.iniEquilibrium(superGeometry, 2, rho, u); sLattice.defineRhoU(superGeometry, 3, rho, u); sLattice.iniEquilibrium(superGeometry, 3, rho, u); sLattice.initialize(); - clout << "Prepare Lattice ... OK" << std::endl; + clout << "Prepare coarse lattice ... OK" << std::endl; } void prepareFineLattice(UnitConverter const& converter, @@ -139,7 +131,7 @@ void prepareFineLattice(UnitConverter const& converter, SuperGeometry2D& superGeometry) { OstreamManager clout(std::cout,"prepareLattice"); - clout << "Prepare Lattice ..." << std::endl; + clout << "Prepare fine lattice ..." << std::endl; const T omega = converter.getLatticeRelaxationFrequency(); @@ -150,25 +142,28 @@ void prepareFineLattice(UnitConverter const& converter, sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega); - T Lx = converter.getLatticeLength(lx); - T Ly = converter.getLatticeLength(ly); - - T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); + const T Lx = converter.getLatticeLength(lx); + const T Ly = converter.getLatticeLength(ly); + const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1); - AnalyticalConst2D iniRho(1); - AnalyticalConst2D iniU(std::vector {0.,0.}); + const T maxVelocity = converter.getCharLatticeVelocity(); + const T radius = ly/2; + std::vector axisPoint{lx, ly/2}; + std::vector axisDirection{1, 0}; + Poiseuille2D u(axisPoint, axisDirection, maxVelocity, radius); + + sLattice.defineRhoU(superGeometry, 1, rho, u); + sLattice.iniEquilibrium(superGeometry, 1, rho, u); + sLattice.defineRhoU(superGeometry, 2, rho, u); + sLattice.iniEquilibrium(superGeometry, 2, rho, u); - sLattice.defineRhoU(superGeometry, 1, iniRho, iniU); - sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU); - sLattice.defineRhoU(superGeometry, 2, iniRho, iniU); - sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU); - sLattice.defineRhoU(superGeometry, 4, iniRho, iniU); - sLattice.iniEquilibrium(superGeometry, 4, iniRho, iniU); + sLattice.defineRhoU(superGeometry, 4, rho, u); + sLattice.iniEquilibrium(superGeometry, 4, rho, u); sLattice.initialize(); - clout << "Prepare Lattice ... OK" << std::endl; + clout << "Prepare fine lattice ... OK" << std::endl; } void getResults(const std::string& prefix, @@ -214,6 +209,13 @@ private: std::unique_ptr> _lattice; public: + static std::unique_ptr> make(IndicatorF2D& domainF, int resolution, T tau, int re) + { + return std::unique_ptr>( + new Grid(domainF, resolution, tau, re) + ); + } + Grid(IndicatorF2D& domainF, int resolution, T tau, int re): _converter(new UnitConverterFromResolutionAndRelaxationTime( resolution, // resolution: number of voxels per charPhysL @@ -221,7 +223,8 @@ public: T{1}, // charPhysLength: reference length of simulation geometry T{1}, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ T{1./re}, // physViscosity: physical kinematic viscosity in __m^2 / s__ - T{1})), // physDensity: physical density in __kg / m^3__ + T{1}, // physDensity: physical density in __kg / m^3__ + T{1})), _cuboids(new CuboidGeometry2D( domainF, _converter->getConversionFactorLength(), @@ -272,6 +275,17 @@ public: { return 1./getScalingFactor(); } + + std::unique_ptr> refine(IndicatorF2D& domainF) + { + return std::unique_ptr>( + new Grid( + domainF, + 2*getConverter().getResolution(), + 2.0*getConverter().getLatticeRelaxationTime() - 0.5, + getConverter().getReynoldsNumber() + )); + } }; template class DESCRIPTOR> @@ -281,8 +295,8 @@ void computeFeq(const Cell& cell, T fEq[DESCRIPTOR::q]) T u[2] {}; cell.computeRhoU(rho, u); const T uSqr = u[0]*u[0] + u[1]*u[1]; - for (int i=0; i < DESCRIPTOR::q; ++i) { - fEq[i] = lbHelpers::equilibrium(i, rho, u, uSqr); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEq[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr); } } @@ -331,8 +345,8 @@ void storeC2F( T fNeq[DESCRIPTOR::q] {}; computeFneq(coarseLattice.get(x,y), fNeq); - for (int i=0; i < DESCRIPTOR::q; ++i) { - c2f_fneq[y][i] = fNeq[i]; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + c2f_fneq[y][iPop] = fNeq[iPop]; } } } @@ -358,8 +372,8 @@ void interpolateStoredC2F( T fNeq[DESCRIPTOR::q] {}; computeFneq(coarseLattice.get(x,y), fNeq); - for (int i=0; i < DESCRIPTOR::q; ++i) { - c2f_fneq[y][i] = order2interpolation(fNeq[i], c2f_fneq[y][i]); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], c2f_fneq[y][iPop]); } } } @@ -382,12 +396,12 @@ void coupleC2F( T fEq[DESCRIPTOR::q] {}; computeFeq(coarseLattice.get(x,y), fEq); - for (int i=0; i < DESCRIPTOR::q; ++i) { - fineLattice.get(0,2*y)[i] = fEq[i] + coarse.getScalingFactor() * c2f_fneq[y][i]; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineLattice.get(0,2*y)[iPop] = fEq[iPop] + coarse.getScalingFactor() * c2f_fneq[y][iPop]; } } - for (int y=2; y < coarseLattice.getNy()-2; ++y) { + for (int y=1; y < coarseLattice.getNy()-2; ++y) { const T rho = order4interpolation( c2f_rho[y-1], c2f_rho[y], @@ -409,15 +423,15 @@ void coupleC2F( ); const T uSqr = u[0]*u[0] + u[1]*u[1]; - for (int i=0; i < DESCRIPTOR::q; ++i) { - const T fneqi = order4interpolation( - c2f_fneq[y-1][i], - c2f_fneq[y][i], - c2f_fneq[y+1][i], - c2f_fneq[y+2][i] - ); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + const T fneq = order4interpolation( + c2f_fneq[y-1][iPop], + c2f_fneq[y][iPop], + c2f_fneq[y+1][iPop], + c2f_fneq[y+2][iPop] + ); - fineLattice.get(0,1+2*y)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; + fineLattice.get(0,1+2*y)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq; } } @@ -440,56 +454,62 @@ void coupleC2F( ); const T uSqr = u[0]*u[0] + u[1]*u[1]; - for (int i=0; i < DESCRIPTOR::q; ++i) { - const T fneqi = order3interpolation( - c2f_fneq[0][i], - c2f_fneq[1][i], - c2f_fneq[2][i] - ); - fineLattice.get(0,1)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + const T fneq = order3interpolation( + c2f_fneq[0][iPop], + c2f_fneq[1][iPop], + c2f_fneq[2][iPop] + ); + fineLattice.get(0,1)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq; } } { + const int maxY = coarseLattice.getNy()-1; const T rho = order3interpolation( - c2f_rho[coarseLattice.getNy()-1], - c2f_rho[coarseLattice.getNy()-2], - c2f_rho[coarseLattice.getNy()-3] + c2f_rho[maxY ], + c2f_rho[maxY-1], + c2f_rho[maxY-2] ); T u[2] {}; u[0] = order3interpolation( - c2f_u[coarseLattice.getNy()-1][0], - c2f_u[coarseLattice.getNy()-2][0], - c2f_u[coarseLattice.getNy()-3][0] + c2f_u[maxY ][0], + c2f_u[maxY-1][0], + c2f_u[maxY-2][0] ); u[1] = order3interpolation( - c2f_u[coarseLattice.getNy()-1][1], - c2f_u[coarseLattice.getNy()-2][1], - c2f_u[coarseLattice.getNy()-3][1] + c2f_u[maxY ][1], + c2f_u[maxY-1][1], + c2f_u[maxY-2][1] ); const T uSqr = u[0]*u[0] + u[1]*u[1]; - for (int i=0; i < DESCRIPTOR::q; ++i) { - const T fneqi = order3interpolation( - c2f_fneq[coarseLattice.getNy()-1][i], - c2f_fneq[coarseLattice.getNy()-2][i], - c2f_fneq[coarseLattice.getNy()-3][i] - ); - fineLattice.get(0,fineLattice.getNy()-2)[i] = lbHelpers::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + const T fneq = order3interpolation( + c2f_fneq[maxY ][iPop], + c2f_fneq[maxY-1][iPop], + c2f_fneq[maxY-2][iPop] + ); + fineLattice.get(0,fineLattice.getNy()-2)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq; } } } template class DESCRIPTOR> -T computeRestrictedFneq(const BlockLatticeView2D& lattice, int x, int y) +void computeRestrictedFneq(const BlockLatticeView2D& lattice, int x, int y, T restrictedFneq[DESCRIPTOR::q]) { - T sum = T{0}; - for (int i=0; i < DESCRIPTOR::q; ++i) { + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { T fNeq[DESCRIPTOR::q] {}; - computeFneq(lattice.get(x + DESCRIPTOR::c[i][0], y + DESCRIPTOR::c[i][1]), fNeq); - sum += fNeq[i]; + computeFneq(lattice.get(x + DESCRIPTOR::c[iPop][0], y + DESCRIPTOR::c[iPop][1]), fNeq); + + for (int jPop=0; jPop < DESCRIPTOR::q; ++jPop) { + restrictedFneq[jPop] += fNeq[jPop]; + } + } + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + restrictedFneq[iPop] /= DESCRIPTOR::q; } - return sum / DESCRIPTOR::q; } template class DESCRIPTOR> @@ -503,10 +523,11 @@ void coupleF2C(Grid& coarse, Grid& fine) T fEq[DESCRIPTOR::q] {}; computeFeq(fineLattice.get(2,2*y), fEq); - const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y); + T fNeq[DESCRIPTOR::q] {}; + computeRestrictedFneq(fineLattice, 2, 2*y, fNeq); - for (int i=0; i < DESCRIPTOR::q; ++i) { - coarseLattice.get(x,y)[i] = fEq[i] + coarse.getInvScalingFactor() * restrictedFneq; + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + coarseLattice.get(x,y)[iPop] = fEq[iPop] + coarse.getInvScalingFactor() * fNeq[iPop]; } } } @@ -526,101 +547,93 @@ int main(int argc, char* argv[]) Vector fineExtend { 0.5 * lx + coarseDeltaX, ly }; IndicatorCuboid2D fineCuboid(fineExtend, fineOrigin); - Grid coarseGrid( - coarseCuboid, - N, - 0.8, - Re); - Grid fineGrid( - fineCuboid, - 2*coarseGrid.getConverter().getResolution(), - 2.0*coarseGrid.getConverter().getLatticeRelaxationTime() - 0.5, - Re); + auto coarseGrid = Grid::make(coarseCuboid, N, 0.8, Re); + auto fineGrid = coarseGrid->refine(fineCuboid); - prepareGeometry(coarseGrid.getConverter(), coarseGrid.getSuperGeometry()); - prepareGeometry(fineGrid.getConverter(), fineGrid.getSuperGeometry()); + prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry()); + prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry()); Dynamics* coarseBulkDynamics; coarseBulkDynamics = new BGKdynamics( - coarseGrid.getConverter().getLatticeRelaxationFrequency(), + coarseGrid->getConverter().getLatticeRelaxationFrequency(), instances::getBulkMomenta()); - sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice()); + sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid->getSuperLattice()); createLocalBoundaryCondition2D(coarseSBoundaryCondition); prepareCoarseLattice( - coarseGrid.getConverter(), - coarseGrid.getSuperLattice(), + coarseGrid->getConverter(), + coarseGrid->getSuperLattice(), *coarseBulkDynamics, coarseSBoundaryCondition, - coarseGrid.getSuperGeometry()); + coarseGrid->getSuperGeometry()); Dynamics* fineBulkDynamics; fineBulkDynamics = new BGKdynamics( - fineGrid.getConverter().getLatticeRelaxationFrequency(), + fineGrid->getConverter().getLatticeRelaxationFrequency(), instances::getBulkMomenta()); - sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice()); + sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid->getSuperLattice()); createLocalBoundaryCondition2D(fineSBoundaryCondition); prepareFineLattice( - fineGrid.getConverter(), - fineGrid.getSuperLattice(), + fineGrid->getConverter(), + fineGrid->getSuperLattice(), *fineBulkDynamics, fineSBoundaryCondition, - fineGrid.getSuperGeometry()); + fineGrid->getSuperGeometry()); clout << "starting simulation..." << endl; Timer timer( - coarseGrid.getConverter().getLatticeTime(maxPhysT), - coarseGrid.getSuperGeometry().getStatistics().getNvoxel()); + coarseGrid->getConverter().getLatticeTime(maxPhysT), + coarseGrid->getSuperGeometry().getStatistics().getNvoxel()); util::ValueTracer converge( - fineGrid.getConverter().getLatticeTime(physInterval), + fineGrid->getConverter().getLatticeTime(physInterval), residuum); timer.start(); - const auto sizeC2F = coarseGrid.getSuperLattice().getBlockLattice(0).getNy(); + const auto sizeC2F = coarseGrid->getSuperLattice().getBlockLattice(0).getNy(); T c2f_rho[sizeC2F] {}; T c2f_u[sizeC2F][2] {}; T c2f_fneq[sizeC2F][DESCRIPTOR::q] {}; - for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) { + for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) { if (converge.hasConverged()) { clout << "Simulation converged." << endl; break; } - storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); - coarseGrid.getSuperLattice().collideAndStream(); - interpolateStoredC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); + storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq); + coarseGrid->getSuperLattice().collideAndStream(); - fineGrid.getSuperLattice().collideAndStream(); - coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq); + interpolateStoredC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq); + fineGrid->getSuperLattice().collideAndStream(); + coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq); - fineGrid.getSuperLattice().collideAndStream(); - storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq); - coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq); + storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq); + fineGrid->getSuperLattice().collideAndStream(); + coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq); - coupleF2C(coarseGrid, fineGrid); + coupleF2C(*coarseGrid, *fineGrid); getResults( "coarse_", - coarseGrid.getSuperLattice(), - coarseGrid.getConverter(), + coarseGrid->getSuperLattice(), + coarseGrid->getConverter(), iT, - coarseGrid.getSuperGeometry(), + coarseGrid->getSuperGeometry(), timer, converge.hasConverged()); getResults( "fine_", - fineGrid.getSuperLattice(), - fineGrid.getConverter(), + fineGrid->getSuperLattice(), + fineGrid->getConverter(), iT, - fineGrid.getSuperGeometry(), + fineGrid->getSuperGeometry(), timer, converge.hasConverged()); - converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true); + converge.takeValue(fineGrid->getSuperLattice().getStatistics().getAverageEnergy(), true); } timer.stop(); -- cgit v1.2.3