From 5cb634a0c8bf7549322168225bbe65b478a4ad8f Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sat, 29 Dec 2018 12:25:34 +0100 Subject: Setup basic coarse and fine lattices Coupling overlap of one coarse grid width. Coarse grid points intersect fine grid points in this area. --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 373 +++++++++++++++++++++--------- 1 file changed, 258 insertions(+), 115 deletions(-) diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index 2f582ab..9415359 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -43,7 +43,7 @@ typedef double T; const T lx = 2.; // length of the channel const T ly = 1.; // height of the channel -int N = 50; // resolution of the model +int N = 20; // resolution of the model const T Re = 10.; // Reynolds number const T maxPhysT = 20.; // max. simulation time in s, SI unit const T physInterval = 0.25; // interval for the convergence check in s @@ -51,15 +51,16 @@ const T residuum = 1e-5; // residuum for the convergence check const T tuner = 0.99; // for partialSlip only: 0->bounceBack, 1->freeSlip -void prepareGeometry( UnitConverter const& converter, - SuperGeometry2D& superGeometry ) { +void prepareGeometry(UnitConverter const& converter, + SuperGeometry2D& superGeometry) +{ - OstreamManager clout( std::cout,"prepareGeometry" ); + OstreamManager clout(std::cout,"prepareGeometry"); clout << "Prepare Geometry ..." << std::endl; - superGeometry.rename( 0,2 ); + superGeometry.rename(0,2); - superGeometry.rename( 2,1,1,1 ); + superGeometry.rename(2,1,0,1); Vector extend; Vector origin; @@ -69,13 +70,13 @@ void prepareGeometry( UnitConverter const& converter, extend[1] = ly; extend[0] = physSpacing / 2; origin[0] -= physSpacing / 4; - IndicatorCuboid2D inflow( extend, origin ); - superGeometry.rename( 2,3,1,inflow ); + IndicatorCuboid2D inflow(extend, origin); + superGeometry.rename(1,3,1,inflow); // Set material number for outflow origin[0] = lx - physSpacing / 4; - IndicatorCuboid2D outflow( extend, origin ); - superGeometry.rename( 2,4,1,outflow ); + IndicatorCuboid2D outflow(extend, origin); + superGeometry.rename(1,4,1,outflow); // Removes all not needed boundary voxels outside the surface superGeometry.clean(); @@ -88,154 +89,296 @@ void prepareGeometry( UnitConverter const& converter, clout << "Prepare Geometry ... OK" << std::endl; } -void prepareLattice( UnitConverter const& converter, - SuperLattice2D& sLattice, - Dynamics& bulkDynamics, - sOnLatticeBoundaryCondition2D& sBoundaryCondition, - SuperGeometry2D& superGeometry ) { - - OstreamManager clout( std::cout,"prepareLattice" ); - clout << "Prepare Lattice ..." << std::endl; +void prepareCoarseLattice(UnitConverter const& converter, + SuperLattice2D& sLattice, + Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition2D& sBoundaryCondition, + SuperGeometry2D& superGeometry) +{ + OstreamManager clout(std::cout,"prepareLattice"); + clout << "Prepare coarse lattice ..." << std::endl; const T omega = converter.getLatticeRelaxationFrequency(); - sLattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics() ); - sLattice.defineDynamics( superGeometry, 1, &bulkDynamics ); - sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack() ); - sLattice.defineDynamics( superGeometry, 3, &bulkDynamics ); - sLattice.defineDynamics( superGeometry, 4, &bulkDynamics ); + sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk + sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack()); + sLattice.defineDynamics(superGeometry, 3, &bulkDynamics); // inflow + + sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega); + + T Lx = converter.getLatticeLength(lx); + T Ly = converter.getLatticeLength(ly); - sBoundaryCondition.addVelocityBoundary( superGeometry, 3, omega ); - sBoundaryCondition.addPressureBoundary( superGeometry, 4, omega ); + T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); + AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1); - T Lx = converter.getLatticeLength( lx ); - T Ly = converter.getLatticeLength( ly ); + AnalyticalConst2D iniRho(1); + AnalyticalConst2D iniU(std::vector {0.,0.}); - T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/( Ly*Ly ); - AnalyticalLinear2D rho( -p0/lx*DESCRIPTOR::invCs2 , 0 , p0*DESCRIPTOR::invCs2+1 ); + sLattice.defineRhoU(superGeometry, 1, iniRho, iniU); + sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU); + sLattice.defineRhoU(superGeometry, 2, iniRho, iniU); + sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU); T maxVelocity = converter.getCharLatticeVelocity(); T distance2Wall = converter.getConversionFactorLength(); - Poiseuille2D u( superGeometry, 3, maxVelocity, distance2Wall ); - - // Initialize all values of distribution functions to their local equilibrium - 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.defineRhoU( superGeometry, 4, rho, u ); - sLattice.iniEquilibrium( superGeometry, 4, rho, u ); + Poiseuille2D u(superGeometry, 3, maxVelocity, distance2Wall); + + sLattice.defineRhoU(superGeometry, 3, rho, u); + sLattice.iniEquilibrium(superGeometry, 3, rho, u); sLattice.initialize(); clout << "Prepare Lattice ... OK" << std::endl; } -void getResults( SuperLattice2D& sLattice, Dynamics& bulkDynamics, - UnitConverter const& converter, int iT, - SuperGeometry2D& superGeometry, Timer& timer, bool hasConverged ) { +void prepareFineLattice(UnitConverter const& converter, + SuperLattice2D& sLattice, + Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition2D& sBoundaryCondition, + SuperGeometry2D& superGeometry) +{ + OstreamManager clout(std::cout,"prepareLattice"); + clout << "Prepare Lattice ..." << std::endl; - OstreamManager clout( std::cout,"getResults" ); + const T omega = converter.getLatticeRelaxationFrequency(); - SuperVTMwriter2D vtmWriter( "poiseuille2d" ); - SuperLatticePhysVelocity2D velocity( sLattice, converter ); - SuperLatticePhysPressure2D pressure( sLattice, converter ); - vtmWriter.addFunctor( velocity ); - vtmWriter.addFunctor( pressure ); + sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk + sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack()); + sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); // outflow - const int vtmIter = converter.getLatticeTime( maxPhysT/20. ); - const int statIter = converter.getLatticeTime( maxPhysT/20. ); + sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega); - if ( iT==0 ) { - SuperLatticeGeometry2D geometry( sLattice, superGeometry ); - SuperLatticeCuboid2D cuboid( sLattice ); - SuperLatticeRank2D rank( sLattice ); - superGeometry.rename( 0,2 ); - vtmWriter.write( geometry ); - vtmWriter.write( cuboid ); - vtmWriter.write( rank ); + T Lx = converter.getLatticeLength(lx); + T Ly = converter.getLatticeLength(ly); + + 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); + sLattice.defineRhoU(superGeometry, 4, iniRho, iniU); + sLattice.iniEquilibrium(superGeometry, 4, iniRho, iniU); + + sLattice.initialize(); + + clout << "Prepare Lattice ... OK" << std::endl; +} +void getResults(const std::string& prefix, + SuperLattice2D& sLattice, + UnitConverter const& converter, int iT, + SuperGeometry2D& superGeometry, Timer& timer, bool hasConverged) +{ + OstreamManager clout(std::cout,"getResults"); + + SuperVTMwriter2D vtmWriter(prefix + "poiseuille2d"); + SuperLatticePhysVelocity2D velocity(sLattice, converter); + SuperLatticePhysPressure2D pressure(sLattice, converter); + SuperLatticeGeometry2D geometry( sLattice, superGeometry ); + vtmWriter.addFunctor(geometry); + vtmWriter.addFunctor(velocity); + vtmWriter.addFunctor(pressure); + + const int vtmIter = converter.getLatticeTime(maxPhysT/100.); + const int statIter = converter.getLatticeTime(maxPhysT/10.); + + if ( iT==0 ) { vtmWriter.createMasterFile(); } - if ( iT%vtmIter==0 || hasConverged ) { - vtmWriter.write( iT ); - - SuperEuklidNorm2D normVel( velocity ); - BlockReduction2D2D planeReduction( normVel, 600, BlockDataSyncMode::ReduceOnly ); - // write output as JPEG - heatmap::write(planeReduction, iT); + if (iT%vtmIter==0 || hasConverged) { + vtmWriter.write(iT); } - if ( iT%statIter==0 || hasConverged ) { - timer.update( iT ); + if (iT%statIter==0 || hasConverged) { + timer.update(iT); timer.printStep(); - sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) ); + sLattice.getStatistics().print(iT,converter.getPhysTime(iT)); } } -int main( int argc, char* argv[] ) { - olbInit( &argc, &argv ); - singleton::directories().setOutputDir( "./tmp/" ); - OstreamManager clout( std::cout,"main" ); - - UnitConverterFromResolutionAndRelaxationTime const converter( - int {N}, // resolution: number of voxels per charPhysL - (T) 0.8, // latticeRelaxationTime: relaxation time, have to be greater than 0.5! - (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.0 // physDensity: physical density in __kg / m^3__ - ); - converter.print(); - converter.write("poiseuille2d"); - - Vector extend( lx, ly ); - Vector origin; - IndicatorCuboid2D cuboid( extend, origin ); - -#ifdef PARALLEL_MODE_MPI - const int noOfCuboids = singleton::mpi().getSize(); -#else - const int noOfCuboids = 1; -#endif - CuboidGeometry2D cuboidGeometry( cuboid, converter.getConversionFactorLength(), noOfCuboids ); - - HeuristicLoadBalancer loadBalancer( cuboidGeometry ); - SuperGeometry2D superGeometry( cuboidGeometry, loadBalancer, 2 ); +template class DESCRIPTOR> +class Grid { +private: + std::unique_ptr> _converter; + std::unique_ptr> _cuboids; + std::unique_ptr> _balancer; + std::unique_ptr> _geometry; + std::unique_ptr> _lattice; + +public: + Grid(IndicatorF2D& domainF, int resolution, T tau, int re): + _converter(new UnitConverterFromResolutionAndRelaxationTime( + resolution, // resolution: number of voxels per charPhysL + tau, // latticeRelaxationTime: relaxation time, have to be greater than 0.5! + 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__ + _cuboids(new CuboidGeometry2D( + domainF, + _converter->getConversionFactorLength(), + 1)), + _balancer(new HeuristicLoadBalancer( + *_cuboids)), + _geometry(new SuperGeometry2D( + *_cuboids, + *_balancer, + 2)), + _lattice(new SuperLattice2D( + *_geometry)) + { + _converter->print(); + } - prepareGeometry( converter, superGeometry ); + UnitConverter& getConverter() + { + return *_converter; + } - SuperLattice2D sLattice( superGeometry ); + CuboidGeometry2D& getCuboidGeometry() + { + return *_cuboids; + } - Dynamics* bulkDynamics; - bulkDynamics = new BGKdynamics( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta() ); + LoadBalancer& getLoadBalancer() + { + return *_balancer; + } - sOnLatticeBoundaryCondition2D sBoundaryCondition( sLattice ); - createInterpBoundaryCondition2D ( sBoundaryCondition ); + SuperGeometry2D& getSuperGeometry() + { + return *_geometry; + } - prepareLattice( converter, sLattice, *bulkDynamics, sBoundaryCondition, superGeometry ); + SuperLattice2D& getSuperLattice() + { + return *_lattice; + } +}; + +int main(int argc, char* argv[]) +{ + olbInit(&argc, &argv); + singleton::directories().setOutputDir("./tmp/"); + OstreamManager clout(std::cout,"main"); + + const T coarseDeltaX = 1./N; + + Vector coarseOrigin { 0.0 , 0.0 }; + Vector coarseExtend { 0.5 * lx, ly }; + IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin); + Vector fineOrigin { 0.5 * lx - coarseDeltaX, 0.0 }; + 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); + + prepareGeometry(coarseGrid.getConverter(), coarseGrid.getSuperGeometry()); + + prepareGeometry(fineGrid.getConverter(), fineGrid.getSuperGeometry()); + + Dynamics* coarseBulkDynamics; + coarseBulkDynamics = new BGKdynamics( + coarseGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice()); + createLocalBoundaryCondition2D(coarseSBoundaryCondition); + + prepareCoarseLattice( + coarseGrid.getConverter(), + coarseGrid.getSuperLattice(), + *coarseBulkDynamics, + coarseSBoundaryCondition, + coarseGrid.getSuperGeometry()); + + Dynamics* fineBulkDynamics; + fineBulkDynamics = new BGKdynamics( + fineGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice()); + createLocalBoundaryCondition2D(fineSBoundaryCondition); + + prepareFineLattice( + fineGrid.getConverter(), + fineGrid.getSuperLattice(), + *fineBulkDynamics, + fineSBoundaryCondition, + fineGrid.getSuperGeometry()); clout << "starting simulation..." << endl; - Timer timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() ); - util::ValueTracer converge( converter.getLatticeTime( physInterval ), residuum ); + Timer timer( + coarseGrid.getConverter().getLatticeTime(maxPhysT), + coarseGrid.getSuperGeometry().getStatistics().getNvoxel()); + util::ValueTracer converge( + coarseGrid.getConverter().getLatticeTime(physInterval), + residuum); timer.start(); - for ( int iT = 0; iT < converter.getLatticeTime( maxPhysT ); ++iT ) { - if ( converge.hasConverged() ) { + for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); iT += 2) { + if (converge.hasConverged()) { clout << "Simulation converged." << endl; - getResults( sLattice, *bulkDynamics, converter, iT, superGeometry, timer, converge.hasConverged() ); - + getResults( + "coarse_", + coarseGrid.getSuperLattice(), + coarseGrid.getConverter(), + iT, + coarseGrid.getSuperGeometry(), + timer, + converge.hasConverged()); break; } - sLattice.collideAndStream(); - - getResults( sLattice, *bulkDynamics, converter, iT, superGeometry, timer, converge.hasConverged() ); - converge.takeValue( sLattice.getStatistics().getAverageEnergy(), true ); + coarseGrid.getSuperLattice().collideAndStream(); + getResults( + "coarse_", + coarseGrid.getSuperLattice(), + coarseGrid.getConverter(), + iT / 2, + coarseGrid.getSuperGeometry(), + timer, + converge.hasConverged()); + + fineGrid.getSuperLattice().collideAndStream(); + getResults( + "fine_", + fineGrid.getSuperLattice(), + fineGrid.getConverter(), + iT, + fineGrid.getSuperGeometry(), + timer, + converge.hasConverged()); + fineGrid.getSuperLattice().collideAndStream(); + getResults( + "fine_", + fineGrid.getSuperLattice(), + fineGrid.getConverter(), + iT+1, + fineGrid.getSuperGeometry(), + timer, + converge.hasConverged()); + + converge.takeValue(coarseGrid.getSuperLattice().getStatistics().getAverageEnergy(), true); } timer.stop(); -- cgit v1.2.3