Setup basic coarse and fine lattices

Coupling overlap of one coarse grid width.
Coarse grid points intersect fine grid points in this area.

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();