C++ MultiBlockLattice2D::periodicity方法代码示例

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

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
           nx, ny, new BGKdynamics<T,DESCRIPTOR>(omega) );

    lattice.periodicity().toggleAll(true); // Set periodic boundaries.

    defineInitialDensityAtCenter(lattice);

    // First part: loop over time iterations.
    for (plint iT=0; iT<maxIter; ++iT) {
        lattice.collideAndStream();
    }

    // Second part: Data analysis.
    Array<T,2> velocity;
    lattice.get(nx/2, ny/2).computeVelocity(velocity);
    pcout << "Velocity in the middle of the lattice: ("
          << velocity[0] << "," << velocity[1] << ")" << endl;

    pcout << "Velocity norm along a horizontal line: " << endl;
    Box2D line(0, 100, ny/2, ny/2);
    pcout << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

    plb_ofstream ofile("profile.dat");
    ofile << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

    pcout << "Average density in the domain: " << computeAverageDensity(lattice) << endl;
    pcout << "Average energy in the domain: " << computeAverageEnergy(lattice) << endl;
}

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

    // For the choice of the parameters G, rho0, and psi0, we refer to the book
    //   Michael C. Sukop and Daniel T. Thorne (2006), 
    //   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
    //   Springer-Verlag Berlin/Heidelberg.
    
    const T omega = 1.0;
    const int nx   = 400;
    const int ny   = 400;
    const T G      = -120.0;
    const int maxIter  = 100001;
    const int saveIter = 100;
    const int statIter = 100;
    
    const T rho0 = 200.0;
    const T deltaRho = 1.0;
    const T psi0 = 4.0;

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
            nx,ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega) );
            
    lattice.periodicity().toggleAll(true);

    // Use a random initial condition, to activate the phase separation.
    applyProcessingFunctional(new RandomInitializer<T,DESCRIPTOR>(rho0,deltaRho), 
                              lattice.getBoundingBox(),lattice);

    // Add the data processor which implements the Shan/Chen interaction potential.
    plint processorLevel = 1;
    integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G, new interparticlePotential::PsiShanChen94<T>(psi0,rho0) ),
            lattice.getBoundingBox(), lattice, processorLevel );

    lattice.initialize();
    
    pcout << "Starting simulation" << endl;
    for (int iT=0; iT<maxIter; ++iT) {
        if (iT%statIter==0) {
            auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
            pcout << iT << ": Average rho fluid one = " << computeAverage(*rho) << endl;
            pcout << "Minimum density: " << computeMin(*rho) << endl;
            pcout << "Maximum density: " << computeMax(*rho) << endl;
        }
        if (iT%saveIter == 0) {
            ImageWriter<T>("leeloo").writeScaledGif (
                    createFileName("rho", iT, 6), *computeDensity(lattice) );
        }

        lattice.collideAndStream();
    }
}

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

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

    IncomprFlowParam<T> parameters(
            (T) 1e-2,  // uMax
            (T) 300.,  // Re
            100,       // N
            5.,        // lx
            1.         // ly 
    );
    const T logT     = (T)0.02;
    const T imSave   = (T)0.1;
    const T vtkSave  = (T)3.;
    const T maxT     = (T)10.1;

    writeLogFile(parameters, "Poiseuille flow");

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

    lattice.periodicity().toggle(0, false);

    OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
        //boundaryCondition = createInterpBoundaryCondition2D<T,DESCRIPTOR>();
        boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

    defineCylinderGeometry(lattice, parameters);
    setupInletAndBulk(lattice, parameters, *boundaryCondition);
    //copyUnknownOnOutlet(lattice, parameters, *boundaryCondition);
    velocityNeumannOutlet(lattice, parameters, *boundaryCondition);
    lattice.initialize();

    // Main loop over time iterations.
    for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
        if ((iT+1)%parameters.nStep(logT)==0) {
            pcout << computeAverageDensity(lattice) << endl;
            pcout << computeAverageEnergy(lattice) << endl;
        }
        if (iT%parameters.nStep(logT)==0) {
            pcout << "step " << iT
                  << "; lattice time=" << lattice.getTimeCounter().getTime()
                  << "; t=" << iT*parameters.getDeltaT()
                  << "; av energy="
                  << setprecision(10) << getStoredAverageEnergy<T>(lattice)
                  << "; av rho="
                  << getStoredAverageDensity<T>(lattice) << endl;
        }

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

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

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

    delete boundaryCondition;
}

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

    if (argc != 2)
    {
        pcout << "Error : Wrong parameters specified." << endl;
        pcout << "1 : N." << endl;
        exit(1);
    }

    const plint N = atoi(argv[1]);

    const T Re = 1.0;
    const T alpha = 1.0; // womersley number

    const plint Nref = 10;

    const T uMaxRef = 0.01;

    const T uMax = uMaxRef /(T)N * (T)Nref; // needed to avoid compressibility errors.

    const T lx  = 100.0;
    const T ly  = 1.0;
    pcout << "uMaxRef=" << uMaxRef << std::endl;
    pcout << "uMax=" << uMax << std::endl;

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

    IncomprFlowParam<T> parameters(uMax, Re, N, lx, ly);

//     The frequency of the force (lattice units)
    T frequency = (T)4*alpha*alpha*parameters.getLatticeNu()
                  / (T)(parameters.getResolution()*parameters.getResolution());

//     The amplitude of the forcing term (lattice units)
    T amplitude = 8. * parameters.getLatticeNu() * parameters.getLatticeU()
                  / ( (T)(parameters.getResolution()*parameters.getResolution()) );

//     Period of the force (lattice units)
    plint tPeriod = (plint)((T)2*pi/frequency + 0.5);

    writeLogFile(parameters,"palabos.log");

    plint nx = parameters.getNx();
    plint ny = parameters.getNy();

    T omega = parameters.getOmega();

    MultiBlockLattice2D<T, NSDESCRIPTOR> lattice (
        nx,ny,new DYNAMICS );

    OnLatticeBoundaryCondition2D<T,NSDESCRIPTOR>*
    boundaryCondition = createLocalBoundaryCondition2D<T,NSDESCRIPTOR>();

    lattice.periodicity().toggle(0,true);

    channelSetup( lattice, parameters, *boundaryCondition, alpha, frequency, amplitude);

    pcout << "Starting simulation" << endl;

    const plint maxIter = tPeriod * 100;
    //const plint tSave = tPeriod / 24;


    T error = T();

    lattice.getTimeCounter().resetTime(1);

    pcout << "Omega = " << omega << ", it period = " << tPeriod << endl;

    util::ValueTracer<T> converge(uMax,N,1.0e-3);
    plint iT = 0;
    for (iT = 0; iT < maxIter; ++iT) {
//         Updating the force in the whole domain
        Array<T,NSDESCRIPTOR<T>::d> force(womersleyForce((T)iT, amplitude, frequency, parameters),0.);
        setExternalVector(lattice,lattice.getBoundingBox(),
                          NSDESCRIPTOR<T>::ExternalField::forceBeginsAt,force);

        T errorTemp = computeRMSerror( lattice,parameters,alpha,iT);
        error += errorTemp;

        //if (iT % tSave == 0) {
        //    pcout << "Writing Gif at time : " << iT << std::endl;
        //    writeGif(lattice,iT);
        //}

        if (iT % tPeriod == 0)
        {
//             The error is averaged over one period
            error /= (T)(tPeriod);
            pcout << "For N = " << N << ", Error = " << error << endl;
            converge.takeValue(error,true);
            if (converge.hasConverged())
            {
                cout << "Simulation converged!\n";
                break;
            }
            error = T();
        }
//.........这里部分代码省略.........