summaryrefslogtreecommitdiff
path: root/apps/adrian/poiseuille2d
diff options
context:
space:
mode:
Diffstat (limited to 'apps/adrian/poiseuille2d')
-rw-r--r--apps/adrian/poiseuille2d/poiseuille2d.cpp373
1 files 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<T,DESCRIPTOR> const& converter,
- SuperGeometry2D<T>& superGeometry ) {
+void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
+ SuperGeometry2D<T>& 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<T,2> extend;
Vector<T,2> origin;
@@ -69,13 +70,13 @@ void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
extend[1] = ly;
extend[0] = physSpacing / 2;
origin[0] -= physSpacing / 4;
- IndicatorCuboid2D<T> inflow( extend, origin );
- superGeometry.rename( 2,3,1,inflow );
+ IndicatorCuboid2D<T> inflow(extend, origin);
+ superGeometry.rename(1,3,1,inflow);
// Set material number for outflow
origin[0] = lx - physSpacing / 4;
- IndicatorCuboid2D<T> outflow( extend, origin );
- superGeometry.rename( 2,4,1,outflow );
+ IndicatorCuboid2D<T> 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<T,DESCRIPTOR> const& converter,
clout << "Prepare Geometry ... OK" << std::endl;
}
-void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
- SuperLattice2D<T, DESCRIPTOR>& sLattice,
- Dynamics<T, DESCRIPTOR>& bulkDynamics,
- sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
- SuperGeometry2D<T>& superGeometry ) {
-
- OstreamManager clout( std::cout,"prepareLattice" );
- clout << "Prepare Lattice ..." << std::endl;
+void prepareCoarseLattice(UnitConverter<T,DESCRIPTOR> const& converter,
+ SuperLattice2D<T, DESCRIPTOR>& sLattice,
+ Dynamics<T, DESCRIPTOR>& bulkDynamics,
+ sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
+ SuperGeometry2D<T>& superGeometry)
+{
+ OstreamManager clout(std::cout,"prepareLattice");
+ clout << "Prepare coarse lattice ..." << std::endl;
const T omega = converter.getLatticeRelaxationFrequency();
- sLattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );
- sLattice.defineDynamics( superGeometry, 1, &bulkDynamics );
- sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>() );
- sLattice.defineDynamics( superGeometry, 3, &bulkDynamics );
- sLattice.defineDynamics( superGeometry, 4, &bulkDynamics );
+ sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
+ sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
+ sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
+ 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<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);
- T Lx = converter.getLatticeLength( lx );
- T Ly = converter.getLatticeLength( ly );
+ AnalyticalConst2D<T,T> iniRho(1);
+ AnalyticalConst2D<T,T> iniU(std::vector<T> {0.,0.});
- T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/( Ly*Ly );
- AnalyticalLinear2D<T,T> rho( -p0/lx*DESCRIPTOR<T>::invCs2 , 0 , p0*DESCRIPTOR<T>::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<T> 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<T> 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<T,DESCRIPTOR>& sLattice, Dynamics<T, DESCRIPTOR>& bulkDynamics,
- UnitConverter<T,DESCRIPTOR> const& converter, int iT,
- SuperGeometry2D<T>& superGeometry, Timer<T>& timer, bool hasConverged ) {
+void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
+ SuperLattice2D<T, DESCRIPTOR>& sLattice,
+ Dynamics<T, DESCRIPTOR>& bulkDynamics,
+ sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
+ SuperGeometry2D<T>& superGeometry)
+{
+ OstreamManager clout(std::cout,"prepareLattice");
+ clout << "Prepare Lattice ..." << std::endl;
- OstreamManager clout( std::cout,"getResults" );
+ const T omega = converter.getLatticeRelaxationFrequency();
- SuperVTMwriter2D<T> vtmWriter( "poiseuille2d" );
- SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity( sLattice, converter );
- SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure( sLattice, converter );
- vtmWriter.addFunctor( velocity );
- vtmWriter.addFunctor( pressure );
+ sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
+ sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
+ sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
+ 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<T, DESCRIPTOR> geometry( sLattice, superGeometry );
- SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice );
- SuperLatticeRank2D<T, DESCRIPTOR> 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<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);
+
+ AnalyticalConst2D<T,T> iniRho(1);
+ AnalyticalConst2D<T,T> iniU(std::vector<T> {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<T,DESCRIPTOR>& sLattice,
+ UnitConverter<T,DESCRIPTOR> const& converter, int iT,
+ SuperGeometry2D<T>& superGeometry, Timer<T>& timer, bool hasConverged)
+{
+ OstreamManager clout(std::cout,"getResults");
+
+ SuperVTMwriter2D<T> vtmWriter(prefix + "poiseuille2d");
+ SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);
+ SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);
+ SuperLatticeGeometry2D<T, DESCRIPTOR> 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<T, DESCRIPTOR> normVel( velocity );
- BlockReduction2D2D<T> 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<T, DESCRIPTOR> 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<T,2> extend( lx, ly );
- Vector<T,2> origin;
- IndicatorCuboid2D<T> cuboid( extend, origin );
-
-#ifdef PARALLEL_MODE_MPI
- const int noOfCuboids = singleton::mpi().getSize();
-#else
- const int noOfCuboids = 1;
-#endif
- CuboidGeometry2D<T> cuboidGeometry( cuboid, converter.getConversionFactorLength(), noOfCuboids );
-
- HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );
- SuperGeometry2D<T> superGeometry( cuboidGeometry, loadBalancer, 2 );
+template <typename T, template<typename> class DESCRIPTOR>
+class Grid {
+private:
+ std::unique_ptr<UnitConverter<T,DESCRIPTOR>> _converter;
+ std::unique_ptr<CuboidGeometry2D<T>> _cuboids;
+ std::unique_ptr<LoadBalancer<T>> _balancer;
+ std::unique_ptr<SuperGeometry2D<T>> _geometry;
+ std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice;
+
+public:
+ Grid(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
+ _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
+ 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<T>(
+ domainF,
+ _converter->getConversionFactorLength(),
+ 1)),
+ _balancer(new HeuristicLoadBalancer<T>(
+ *_cuboids)),
+ _geometry(new SuperGeometry2D<T>(
+ *_cuboids,
+ *_balancer,
+ 2)),
+ _lattice(new SuperLattice2D<T,DESCRIPTOR>(
+ *_geometry))
+ {
+ _converter->print();
+ }
- prepareGeometry( converter, superGeometry );
+ UnitConverter<T,DESCRIPTOR>& getConverter()
+ {
+ return *_converter;
+ }
- SuperLattice2D<T, DESCRIPTOR> sLattice( superGeometry );
+ CuboidGeometry2D<T>& getCuboidGeometry()
+ {
+ return *_cuboids;
+ }
- Dynamics<T, DESCRIPTOR>* bulkDynamics;
- bulkDynamics = new BGKdynamics<T, DESCRIPTOR>( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>() );
+ LoadBalancer<T>& getLoadBalancer()
+ {
+ return *_balancer;
+ }
- sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> sBoundaryCondition( sLattice );
- createInterpBoundaryCondition2D<T, DESCRIPTOR> ( sBoundaryCondition );
+ SuperGeometry2D<T>& getSuperGeometry()
+ {
+ return *_geometry;
+ }
- prepareLattice( converter, sLattice, *bulkDynamics, sBoundaryCondition, superGeometry );
+ SuperLattice2D<T,DESCRIPTOR>& 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<T,2> coarseOrigin { 0.0 , 0.0 };
+ Vector<T,2> coarseExtend { 0.5 * lx, ly };
+ IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
+ Vector<T,2> fineOrigin { 0.5 * lx - coarseDeltaX, 0.0 };
+ Vector<T,2> fineExtend { 0.5 * lx + coarseDeltaX, ly };
+ IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
+
+ Grid<T,DESCRIPTOR> coarseGrid(
+ coarseCuboid,
+ N,
+ 0.8,
+ Re);
+ Grid<T,DESCRIPTOR> 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<T, DESCRIPTOR>* coarseBulkDynamics;
+ coarseBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
+ coarseGrid.getConverter().getLatticeRelaxationFrequency(),
+ instances::getBulkMomenta<T, DESCRIPTOR>());
+
+ sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
+ createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);
+
+ prepareCoarseLattice(
+ coarseGrid.getConverter(),
+ coarseGrid.getSuperLattice(),
+ *coarseBulkDynamics,
+ coarseSBoundaryCondition,
+ coarseGrid.getSuperGeometry());
+
+ Dynamics<T, DESCRIPTOR>* fineBulkDynamics;
+ fineBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
+ fineGrid.getConverter().getLatticeRelaxationFrequency(),
+ instances::getBulkMomenta<T, DESCRIPTOR>());
+
+ sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
+ createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
+
+ prepareFineLattice(
+ fineGrid.getConverter(),
+ fineGrid.getSuperLattice(),
+ *fineBulkDynamics,
+ fineSBoundaryCondition,
+ fineGrid.getSuperGeometry());
clout << "starting simulation..." << endl;
- Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
- util::ValueTracer<T> converge( converter.getLatticeTime( physInterval ), residuum );
+ Timer<T> timer(
+ coarseGrid.getConverter().getLatticeTime(maxPhysT),
+ coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
+ util::ValueTracer<T> 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();