C++ IncomprFlowParam类代码示例

void cavitySetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
                  IncomprFlowParam<T> const& parameters,
                  OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();
    const plint nz = parameters.getNz();
    Box3D topLid = Box3D(0, nx-1, ny-1, ny-1, 0, nz-1);
    Box3D bottomLid = Box3D(0, nx-1, 0, 0, 0, nz-1);
    Box3D everythingButTopLid = Box3D(0, nx-1, 0, ny-2, 0, nz-1);

    /*
    instantiateOuterNLDboundary(lattice, lattice.getBoundingBox());
    setOuterNLDboundaryDynamics(lattice, lattice.getBackgroundDynamics().clone(),
                                lattice.getBoundingBox(), boundary::dirichlet);
    setOuterNLDboundaryDynamics(lattice, lattice.getBackgroundDynamics().clone(), bottomLid, boundary::neumann);
    defineDynamics(lattice, bottomLid, lattice.getBackgroundDynamics().clone());
    */

    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, lattice.getBoundingBox(), boundary::dirichlet);

    T u = sqrt((T)2)/(T)2 * parameters.getLatticeU();
    initializeAtEquilibrium(lattice, everythingButTopLid, (T) 1., Array<T,3>(0.,0.,0.) );
    initializeAtEquilibrium(lattice, topLid, (T) 1., Array<T,3>(u,0.,u) );
    setBoundaryVelocity(lattice, topLid, Array<T,3>(u,0.,u) );

    lattice.initialize();
}

void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
			  IncomprFlowParam<T> const& parameters, plint iter)
{
	T dx = parameters.getDeltaX();
	T dt = parameters.getDeltaT();
	VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 6), dx);
	vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
	vtkOut.writeData<2,float>(*computeVelocity(lattice), "velocity", dx/dt);
}

void writeVTK(BlockLatticeT& lattice,
              IncomprFlowParam<T> const& parameters, plint iter)
{
    T dx = parameters.getDeltaX();
    T dt = parameters.getDeltaT();
    VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);
    vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
    vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);
    vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), "vorticity", 1./dt);
}

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

    boundaryCondition.addVelocityBoundary0P (
            Box2D(nx-1,nx-1, 0,ny-1), lattice, boundary::outflow );
}

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

    // On the right boundary, we copy the unknown populations from previous locations
    integrateProcessingFunctional(new CopyUnknownPopulationsFunctional2D<T,DESCRIPTOR, 0, +1>,
                                  Box2D(nx-1,nx-1, 0,ny-1), lattice);
}

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

    int cx     = nx/4;
    int cy     = ny/2+ny/10;
    int radius = cy/4;

    createCylinder(lattice, cx, cy, radius);
}

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

/// Velocity on the parabolic Womersley profile
T womersleyVelocity(plint iY, T t, T A, T omega, T alpha, IncomprFlowParam<T> const& parameters)
{
    const complex<T> I(0.,1.);
    T y = (T)iY / parameters.getResolution();

    return ( A / (I * omega) * std::exp(I * omega * t) *
             ( (T)1. - std::cosh(std::sqrt((T)2.)*(y-(T).5)*(alpha+I*alpha)) /
               std::cosh(std::sqrt((T)2.)/(T)2. * (alpha+I*alpha)) ) ).real();
}

void writeGifs(BlockLatticeT& lattice,
               IncomprFlowParam<T> const& parameters, plint iter)
{
    const plint imSize = 600;
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();
    const plint nz = parameters.getNz();
    const plint zComponent = 2;

    Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);
    ImageWriter<T> imageWriter("leeloo");

    imageWriter.writeScaledGif( createFileName("uz", iter, 6),
                                *computeVelocityComponent (lattice, slice, zComponent),
                                imSize, imSize );
    imageWriter.writeScaledGif( createFileName("uNorm", iter, 6),
                                *computeVelocityNorm (lattice, slice),
                                imSize, imSize );
}

void cavitySetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
                  IncomprFlowParam<T> const& parameters,
                  OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
    const plint nx = parameters.getNx();
    const plint ny = parameters.getNy();
    const plint nz = parameters.getNz();
    Box3D topLid = Box3D(0, nx-1, ny-1, ny-1, 0, nz-1);
    Box3D everythingButTopLid = Box3D(0, nx-1, 0, ny-2, 0, nz-1);

    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    T u = std::sqrt((T)2)/(T)2 * parameters.getLatticeU();
    initializeAtEquilibrium(lattice, everythingButTopLid, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );
    initializeAtEquilibrium(lattice, topLid, (T)1., Array<T,3>(u,(T)0.,u) );
    setBoundaryVelocity(lattice, topLid, Array<T,3>(u,(T)0.,u) );

    lattice.initialize();
}

/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder
void cylinderSetup( BlockLatticeBase2D<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, 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<T>(parameters) );
	setBoundaryDensity (
		lattice, outlet,
		ConstantDensity<T>(1.) );
	initializeAtEquilibrium (
		lattice, lattice.getBoundingBox(),
		PoiseuilleVelocityAndDensity<T>(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<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;
}

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

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

    // Create Velocity boundary conditions on inlet
    boundaryCondition.addVelocityBoundary0N(Box2D(   0,   0,   0,ny-1), lattice);

    setBoundaryVelocity (
            lattice, Box2D(   0,   0,   0,ny-1),
            PoiseuilleVelocity<T>(parameters) );
    initializeAtEquilibrium (
            lattice, lattice.getBoundingBox(),
            PoiseuilleVelocityAndDensity<T,DESCRIPTOR>(parameters) );
}