C++ MultiBlockLattice3D::collideAndStream方法代码示例

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);

    global::directories().setOutputDir("./tmp/");

    // Use the class IncomprFlowParam to convert from
    //   dimensionless variables to lattice units, in the
    //   context of incompressible flows.
    IncomprFlowParam<T> parameters(
        (T) 1e-2,  // Reference velocity (the maximum velocity
        //   in the Poiseuille profile) in lattice units.
        (T) 100.,  // Reynolds number
        30,        // Resolution of the reference length (channel height).
        3.,        // Channel length in dimensionless variables
        1.,        // Channel height in dimensionless variables
        1.         // Channel depth in dimensionless variables
    );
    const T imSave   = (T)0.02; // Time intervals at which to save GIF
    //   images, in dimensionless time units.
    const T maxT     = (T)2.5;  // Total simulation time, in dimensionless
    //   time units.

    writeLogFile(parameters, "3D Poiseuille flow");

    MultiBlockLattice3D<T, DESCRIPTOR> lattice (
        parameters.getNx(), parameters.getNy(), parameters.getNz(),
        new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
    //boundaryCondition = createInterpBoundaryCondition3D<T,DESCRIPTOR>();
    boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

    channelSetup(lattice, parameters, *boundaryCondition);

    // Main loop over time iterations.
    for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
        if (iT%parameters.nStep(imSave)==0) {
            pcout << "Saving Gif at time step " << iT << endl;
            writeGifs(lattice, iT);
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }

    delete boundaryCondition;
}

int main(int argc, char* argv[]) {
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // uMax
            (T) 100.,  // Re
            100,       // N
            1.,        // lx
            1.,        // ly
            1.         // lz
    );
    const plint logIter        = 10;
    const plint imageIter      = 60;
    const plint checkPointIter = 200;
    plint iniT, endT;

    if (argc != 3) {
        pcout << "Error; the syntax is \"" << argv[0] << " start-iter end-iter\"," << endl;
        return -1;
    }

    stringstream iniTstr, endTstr;
    iniTstr << argv[1]; iniTstr >> iniT;
    endTstr << argv[2]; endTstr >> endT;

    writeLogFile(parameters, "3D diagonal cavity");

    MultiBlockLattice3D<T, DESCRIPTOR> lattice (
            parameters.getNx(), parameters.getNy(), parameters.getNz(),
            new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
        //= createInterpBoundaryCondition3D<T,DESCRIPTOR>();
        = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

    cavitySetup(lattice, parameters, *boundaryCondition);
    delete boundaryCondition;

    // Load saved data from a previous simulation, if the initial time step
    //   is larger than zero.
    if (iniT>0) {
        //loadRawMultiBlock(lattice, "checkpoint.dat");
        loadBinaryBlock(lattice, "checkpoint.dat");
    }

    // Main loop over time iterations.
    for (plint iT=iniT; iT<endT; ++iT) {
        if (iT%imageIter==0) {
            pcout << "Writing Gif ..." << endl;
            writeGifs(lattice, parameters, iT);
        }

        if (iT%checkPointIter==0 && iT>iniT) {
            pcout << "Saving the state of the simulation ..." << endl;
            //saveRawMultiBlock(lattice, "checkpoint.dat");
            saveBinaryBlock(lattice, "checkpoint.dat");
        }

        // Lattice Boltzmann iteration step.
        lattice.collideAndStream();

        if (iT%logIter==0) {
            pcout << "step " << iT
                  << "; t=" << iT*parameters.getDeltaT()
                  << "; av energy="
                  << setprecision(10) << getStoredAverageEnergy<T>(lattice)
                  << "; av rho="
                  << setprecision(10) << getStoredAverageDensity<T>(lattice) << endl;
        }
    }
}

int main(int argc, char* argv[]) {

    plbInit(&argc, &argv);
    //defaultMultiBlockPolicy3D().toggleBlockingCommunication(true);

    plint N;
    try {
        global::argv(1).read(N);
    }
    catch(...)
    {
        pcout << "Wrong parameters. The syntax is " << std::endl;
        pcout << argv[0] << " N" << std::endl;
        pcout << "where N is the resolution. The benchmark cases published " << std::endl;
        pcout << "on the Palabos Wiki use N=100, N=400, N=1000, or N=4000." << std::endl;
        exit(1);
    }

    pcout << "Starting benchmark with " << N+1 << "x" << N+1 << "x" << N+1 << " grid points "
          << "(approx. 2 minutes on modern processors)." << std::endl;


    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // uMax
            (T) 1.,    // Re
            N,         // N
            1.,        // lx
            1.,        // ly
            1.         // lz
    );


    MultiBlockLattice3D<T, DESCRIPTOR> lattice (
            parameters.getNx(), parameters.getNy(), parameters.getNz(),
            new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

    plint numCores = global::mpi().getSize();
    pcout << "Number of MPI threads: " << numCores << std::endl;
    // Current cores run approximately at 5 Mega Sus.
    T estimateSus= 5.e6*numCores;
    // The benchmark should run for approximately two minutes
    // (2*60 seconds).
    T wishNumSeconds = 60.;
    plint numCells = lattice.getBoundingBox().nCells();

    // Run at least three iterations.
    plint numIter = std::max( (plint)3,
                              (plint)(estimateSus*wishNumSeconds/numCells+0.5));

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
        = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

    cavitySetup(lattice, parameters, *boundaryCondition);

    // Run the benchmark once "to warm up the machine".
    for (plint iT=0; iT<numIter; ++iT) {
        lattice.collideAndStream();
    }

    // Run the benchmark for good.
    global::timer("benchmark").start();
    global::profiler().turnOn();
    for (plint iT=0; iT<numIter; ++iT) {
        lattice.collideAndStream();
    }

    pcout << "After " << numIter << " iterations: "
          << (T) (numCells*numIter) /
             global::timer("benchmark").getTime() / 1.e6
          << " Mega site updates per second." << std::endl << std::endl;

    global::profiler().writeReport();

    delete boundaryCondition;
}

int main(int argc, char* argv[]) {

    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // uMax
            (T) 10.,   // Re
            30,        // N
            1.,        // lx
            1.,        // ly
            1.         // lz
    );
    const T logT     = (T)1/(T)100;
    const T imSave   = (T)1/(T)10;
    const T vtkSave  = (T)1;
    const T maxT     = (T)10.1;

    pcout << "omega= " << parameters.getOmega() << std::endl;
    writeLogFile(parameters, "3D diagonal cavity");

    T omega = parameters.getOmega();

#ifdef USE_MRT
    plint mrtId = 0;
    mrtParam<T,DESCRIPTOR>().set(mrtId,MRTparam<T,DESCRIPTOR>(omega));
#endif


    MultiBlockLattice3D<T, DESCRIPTOR> lattice (
            parameters.getNx(), parameters.getNy(), parameters.getNz(), DYNAMICS );

#ifdef USE_ASINARI
        integrateProcessingFunctional(new AsinariPostCollide3D<T,DESCRIPTOR>, lattice.getBoundingBox(), lattice, 0);
#endif

    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
        //= createInterpBoundaryCondition3D<T,DESCRIPTOR>();
        = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

    cavitySetup(lattice, parameters, *boundaryCondition);

    T previousIterationTime = T();
    // Loop over main time iteration.
    for (plint iT=0; iT<parameters.nStep(maxT); ++iT) {
        global::timer("mainLoop").restart();

        if (iT%parameters.nStep(imSave)==0) {
            pcout << "Writing Gif ..." << endl;
            writeGifs(lattice, parameters, iT);
        }

        if (iT%parameters.nStep(vtkSave)==0 && iT>0) {
            pcout << "Saving VTK file ..." << endl;
            writeVTK(lattice, parameters, iT);
        }

        if (iT%parameters.nStep(logT)==0) {
            pcout << "step " << iT
                  << "; t=" << iT*parameters.getDeltaT();
        }

        // Execute a time iteration.
        lattice.collideAndStream();

        // Access averages from internal statistics ( their value is defined
        //   only after the call to lattice.collideAndStream() )
        if (iT%parameters.nStep(logT)==0) {
            pcout << "; av energy="
                  << setprecision(10) << getStoredAverageEnergy<T>(lattice)
                  << "; av rho="
                  << setprecision(10) << getStoredAverageDensity<T>(lattice) << endl;
            pcout << "Time spent during previous iteration: "
                  << previousIterationTime << endl;
        }

        previousIterationTime = global::timer("mainLoop").stop();
    }

    delete boundaryCondition;
}

int main(int argc, char* argv[]) {

    /// === 1st Step: Initialization ===

  olbInit(&argc, &argv);
  singleton::directories().setOutputDir("./tmp/");
  OstreamManager clout(std::cout,"main");

  LBconverter<T> converter(
    (int) 3,                               // dim
    (T)   1./20.,                          // latticeL_
    (T)   4e-2,                            // latticeU_
    (T)   1./100.,                         // charNu_
    (T)   1.                               // charL_ = 1,
  );

  writeLogFile(converter, "backwardFacingStep2d");

/// === 3rd Step: Prepare Lattice ===

  //const T maxT     = (T)100.1;
  const T maxT     = (T)0.21;


  writeLogFile(converter, "3D backward facing step");

#ifndef PARALLEL_MODE_MPI  // sequential program execution
  BlockLattice3D<T, DESCRIPTOR> lattice(
    converter.numNodes(18), converter.numNodes(1.5), converter.numNodes(1.5) );
#else                      // parallel program execution
  MultiBlockLattice3D<T, DESCRIPTOR> lattice ( createRegularDataDistribution(
        //converter.numNodes(18), converter.numNodes(1.5), converter.numNodes(1.5) ) );//401*31*31
        converter.numNodes(36), converter.numNodes(12), converter.numNodes(12) ) );
#endif

  ConstRhoBGKdynamics<T, DESCRIPTOR> bulkDynamics (
    converter.getOmega(),
    instances::getBulkMomenta<T,DESCRIPTOR>()
  );

  // choose between local and non-local boundary condition
  OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
  //boundaryCondition = createInterpBoundaryCondition3D(lattice);
  boundaryCondition = createLocalBoundaryCondition3D(lattice);

  prepareLattice(converter, lattice, bulkDynamics, *boundaryCondition);

    /// === 4th Step: Main Loop with Timer ===

	T st, et;
    int iT=0;
      st=time_sec();
//    setBoundaryValues(converter, lattice, iT);
      et=time_sec();
    cout << "set boundary: " << et - st << endl;
      st=time_sec();
    lattice.collideAndStream();
    getResults(lattice, converter, iT, maxT);
  //for (iT=0; iT < converter.numTimeSteps(maxT); ++iT) {
  //for (iT=1; iT < converter.numTimeSteps(maxT); ++iT) {
  for (iT=1; iT < 100; ++iT) {

    /// === 5th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues(converter, lattice, iT);


    /// === 6th Step: Collide and Stream Execution ===
    lattice.collideAndStream();

    /// === 7th Step: Computation and Output of the Results ===
    getResults(lattice, converter, iT, maxT);


  }
      et=time_sec();
    cout << "steps:"<<converter.numTimeSteps(maxT)<<"kernel time: " << et - st << endl;
  std::cout<<"nx:"<<converter.numNodes(36)<<"ny:"<<converter.numNodes(12)<<"nz:"<<converter.numNodes(12)<<std::endl;

  delete boundaryCondition;
}

//.........这里部分代码省略.........
    T rho_s = 1100.;
    T m = r_*r_*r_*4./3.*3.14*rho_s;

    PhysUnits3D<T> units(2.*r_,v_inf,nu_f,lx,ly,lz,N,uMax,rho_f);

    IncomprFlowParam<T> parameters(units.getLbParam());

    plint nx = parameters.getNx(), ny = parameters.getNy(), nz = parameters.getNz();

    // get lattice decomposition from LIGGGHTS and create lattice according to parallelization
    // given in the LIGGGHTS input script
    LatticeDecomposition lDec(parameters.getNx(),parameters.getNy(),parameters.getNz(),
                              wrapper.lmp);
    SparseBlockStructure3D blockStructure = lDec.getBlockDistribution();
    ExplicitThreadAttribution* threadAttribution = lDec.getThreadAttribution();
    plint envelopeWidth = 1;

    MultiBlockLattice3D<T, DESCRIPTOR>
    lattice (MultiBlockManagement3D (blockStructure, threadAttribution, envelopeWidth ),
             defaultMultiBlockPolicy3D().getBlockCommunicator(),
             defaultMultiBlockPolicy3D().getCombinedStatistics(),
             defaultMultiBlockPolicy3D().getMultiCellAccess<T,DESCRIPTOR>(),
             new DYNAMICS );

    defineDynamics(lattice,lattice.getBoundingBox(),new DYNAMICS);


    const T maxT = ceil(3.*lz/v_inf);
    const T vtkT = 0.1;
    const T logT = 0.0000001;

    const plint maxSteps = units.getLbSteps(maxT);
    const plint vtkSteps = max<plint>(units.getLbSteps(vtkT),1);
    const plint logSteps = max<plint>(units.getLbSteps(logT),1);

    writeLogFile(parameters, "sedimenting spheres benchmark");


    lattice.initialize();
    T dt_phys = units.getPhysTime(1);
    plint demSubsteps = 10;
    T dt_dem = dt_phys/(T)demSubsteps;


    pcout << "------------------------------\n"
          << "omega: " << parameters.getOmega() << "\n"
          << "dt_phys: " << dt_phys << "\n"
          << "maxT: " << maxT << " | maxSteps: " << maxSteps << "\n"
          << "v_inf: " << v_inf << "\n"
          << "Re : " << parameters.getRe() << "\n"
          << "vtkT: " << vtkT << " | vtkSteps: " << vtkSteps << "\n"
          << "grid size: " << nx << " " << ny << " " << nz << "\n"
          << "------------------------------" << std::endl;
    // set timestep and output directory
    wrapper.setVariable("t_step",dt_dem);
    wrapper.setVariable("dmp_stp",vtkSteps*demSubsteps);
    wrapper.setVariable("dmp_dir",demOutDir);


    wrapper.execFile("in2.lbdem");
    wrapper.runUpto(demSubsteps-1);

    clock_t start = clock();
    clock_t loop = clock();
    clock_t end = clock();
    // Loop over main time iteration.
    for (plint iT=0; iT<=maxSteps; ++iT) {

        bool initWithVel = false;
        setSpheresOnLattice(lattice,wrapper,units,initWithVel);


        if(iT%vtkSteps == 0 && iT > 0) // LIGGGHTS does not write at timestep 0
            writeVTK(lattice,parameters,units,iT);

        lattice.collideAndStream();

        getForcesFromLattice(lattice,wrapper,units);

        wrapper.run(demSubsteps);


        if(iT%logSteps == 0) {
            end = clock();
            T time = difftime(end,loop)/((T)CLOCKS_PER_SEC);
            T totaltime = difftime(end,start)/((T)CLOCKS_PER_SEC);
            T mlups = ((T) (lattice.getNx()*lattice.getNy()*lattice.getNz()*logSteps))/time/1e6;
            pcout << "time: " << time << " " ;
            pcout << "calculating at " << mlups << " MLU/s"
                  << " | total time running: " << totaltime << std::endl;
            loop = clock();
        }
    }
    T totaltime = difftime(end,start)/((T)CLOCKS_PER_SEC);
    T totalmlups = ((T) (lattice.getNx()*lattice.getNy()*lattice.getNz()*(maxSteps+1)))/totaltime/1e6;
    pcout << " ********************** \n"
          << "total time: " << totaltime
          << " calculating at " << totalmlups << " MLU/s" << std::endl;

}