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

void halfCircleSetup (
        MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint N, plint radius,
        OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
    // The channel is pressure-driven, with a difference deltaRho
    //   between inlet and outlet.
    T deltaRho = 1.e-2;
    T rhoIn  = 1. + deltaRho/2.;
    T rhoOut = 1. - deltaRho/2.;

    Box2D inlet (0,     N/2, N/2, N/2);
    Box2D outlet(N/2+1, N,   N/2, N/2);

    boundaryCondition.addPressureBoundary1P(inlet, lattice);
    boundaryCondition.addPressureBoundary1P(outlet, lattice);

    // Specify the inlet and outlet density.
    setBoundaryDensity (lattice, inlet, rhoIn);
    setBoundaryDensity (lattice, outlet, rhoOut);

    // Create the initial condition.
    Array<T,2> zeroVelocity((T)0.,(T)0.);
    T constantDensity = (T)1;
    initializeAtEquilibrium (
       lattice, lattice.getBoundingBox(), constantDensity, zeroVelocity );

    defineDynamics(lattice, lattice.getBoundingBox(),
                   new BounceBackNodes<T>(N, radius),
                   new BounceBack<T,DESCRIPTOR>);

    lattice.initialize();
}

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
        2.,        // lx
        1.         // ly 
    );

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

    writeLogFile(parameters, "Poiseuille flow");

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
              nx, ny,
              new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );
    OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
        boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();
    createPoiseuilleBoundaries(lattice, parameters, *boundaryCondition);
    lattice.initialize();

    // The following command opens a text-file, in which the velocity-profile
    // in the middle of the channel will be written and several successive
    // time steps. Note the use of plb_ofstream instead of the standard C++
    // ofstream, which is required to guarantee a consistent behavior in MPI-
    // parallel programs.
    plb_ofstream successiveProfiles("velocityProfiles.dat");

    // Main loop over time steps.
    for (plint iT=0; iT<10000; ++iT) {
        if (iT%1000==0) {
            pcout << "At iteration step " << iT
                  << ", the density along the channel is " << endl;
            pcout << setprecision(7)
                  << *computeDensity(lattice, Box2D(0, nx-1, ny/2, ny/2))
                  << endl << endl;

            Box2D profileSection(nx/2, nx/2, 0, ny-1);
            successiveProfiles
                << setprecision(4)
                  // (2) Convert from lattice to physical units.
                << *multiply (
                       parameters.getDeltaX() / parameters.getDeltaT(),
                  // (1) Compute velocity norm along the chosen section.
                       *computeVelocityNorm (lattice, profileSection) )
                << endl;

        }

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

    delete boundaryCondition;
}

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

// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a square sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)
{
    // Initialize constant density everywhere.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), rho0, u0 );

    // And slightly higher density in the central box.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), initializeRhoOnCircle );

    lattice.initialize();
}

/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder
void cylinderSetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
                    IncomprFlowParam<T> const& parameters,
                    OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();

    Box2D outlet(nx-1,nx-1, 1,ny-2);

    // Create Velocity boundary conditions everywhere
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
        lattice, Box2D(0, nx-1, 0, 0) );
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
        lattice, Box2D(0, nx-1, ny-1, ny-1) );
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
        lattice, Box2D(0,0, 1,ny-2) );
    // .. except on right boundary, where we prefer a fixed-pressure condition.
    boundaryCondition.setPressureConditionOnBlockBoundaries (
        lattice, outlet );

    setBoundaryVelocity (
        lattice, lattice.getBoundingBox(),
        PoiseuilleVelocity<T>(parameters) );
    setBoundaryDensity (
        lattice, outlet,
        ConstantDensity<T>(1.) );
    initializeAtEquilibrium (
        lattice, lattice.getBoundingBox(),
        PoiseuilleVelocityAndDensity<T,DESCRIPTOR>(parameters) );

    plint cx = nx/4;
    plint cy = ny/2+2;
    plint r  = cy/4;
    DotList2D cylinderShape;
    for (plint iX=0; iX<nx; ++iX) {
        for (plint iY=0; iY<ny; ++iY) {
            if ( (iX-cx)*(iX-cx) + (iY-cy)*(iY-cy) < r*r ) {
                cylinderShape.addDot(Dot2D(iX,iY));
            }
        }
    }
    defineDynamics(lattice, cylinderShape, new BounceBack<T,DESCRIPTOR>);

    lattice.initialize();
}

void cavitySetup( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
                  IncomprFlowParam<T> const& parameters,
                  OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();

    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), 1., Array<T,2>(0.,0.) );

    T u = parameters.getLatticeU();
    setBoundaryVelocity(lattice, Box2D(1, nx-2, ny-1, ny-1), Array<T,2>(u,0.) );
    initializeAtEquilibrium(lattice, Box2D(1, nx-2, ny-1, ny-1), 1., Array<T,2>(u,0.) );

    lattice.initialize();
}

/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder
void cylinderSetup( MultiBlockLattice2D<PlbT,DESCRIPTOR>& lattice,
                    IncomprFlowParam<PlbT> const& parameters,
                    OnLatticeBoundaryCondition2D<PlbT,DESCRIPTOR>& boundaryCondition )
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();
    Box2D outlet(nx-1,nx-1, 1, ny-2);

    // Create Velocity boundary conditions everywhere
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
            lattice, Box2D(0, 0, 1, ny-2) );
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
            lattice, Box2D(0, nx-1, 0, 0) );
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
            lattice, Box2D(0, nx-1, ny-1, ny-1) );
    // .. except on right boundary, where we prefer an outflow condition
    //    (zero velocity-gradient).
    boundaryCondition.setVelocityConditionOnBlockBoundaries (
            lattice, Box2D(nx-1, nx-1, 1, ny-2), boundary::outflow );

    setBoundaryVelocity (
            lattice, lattice.getBoundingBox(),
            PoiseuilleVelocity<PlbT>(parameters) );
    setBoundaryDensity (
            lattice, outlet,
            ConstantDensity<PlbT>(1.) );
    initializeAtEquilibrium (
            lattice, lattice.getBoundingBox(),
            PoiseuilleVelocityAndDensity<PlbT>(parameters) );

    plint cx     = nx/4;
    plint cy     = ny/2+2; // cy is slightly offset to avoid full symmetry,
                          //   and to get a Von Karman Vortex street.
    plint radius = cy/4;
    defineDynamics(lattice, lattice.getBoundingBox(),
                   new CylinderShapeDomain2D<T>(cx,cy,radius),
                   new plb::BounceBack<PlbT,DESCRIPTOR>);

    lattice.initialize();
}

void channelSetup( MultiBlockLattice2D<T,NSDESCRIPTOR>& lattice,
                   IncomprFlowParam<T> const& parameters,
                   OnLatticeBoundaryCondition2D<T,NSDESCRIPTOR>& boundaryCondition,
                   T alpha, T frequency, T amplitude)
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();

    Box2D bottom(   0,nx-1,   0,   0);
    Box2D top(   0,nx-1,   ny-1,   ny-1);

    boundaryCondition.addVelocityBoundary1N(bottom, lattice);
    boundaryCondition.addPressureBoundary1P(top,    lattice);

    Array<T,2> u((T)0.,(T)0.);
    setBoundaryVelocity( lattice, lattice.getBoundingBox(), u );
    initializeAtEquilibrium(lattice,lattice.getBoundingBox(),(T)1.0,u);

    Array<T,NSDESCRIPTOR<T>::d> force(womersleyForce((T)0, amplitude, frequency, parameters),0.);
    setExternalVector(lattice,lattice.getBoundingBox(),NSDESCRIPTOR<T>::ExternalField::forceBeginsAt,force);

    lattice.initialize();
}

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;
}