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

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 iniLattice( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
                 VoxelizedDomain3D<T>& voxelizedDomain )
{
    // Switch all remaining outer cells to no-dynamics, except the outer
    //   boundary layer, and keep the rest as BGKdynamics.
    defineDynamics(lattice, voxelizedDomain.getVoxelMatrix(), lattice.getBoundingBox(),
                   new NoDynamics<T,DESCRIPTOR>, voxelFlag::outside);
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), (T) 1., Array<T,3>((T)0.,(T)0.,(T)0.));
    lattice.initialize();
}

void channelSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
                   IncomprFlowParam<T> const& parameters,
                   OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
    // Create Velocity boundary conditions
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    setBoundaryVelocity (
        lattice, lattice.getBoundingBox(),
        PoiseuilleVelocity<T>(parameters) );

    lattice.initialize();
}

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

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

    plbInit(&argc, &argv);

    T uMax;

    plint N;

    T nu_f,d_part,v_frac, v_inf;

    std::string outDir;

    try {
        global::argv(1).read(d_part);
        global::argv(2).read(N);
        global::argv(3).read(v_frac);
        global::argv(4).read(nu_f);
        global::argv(5).read(v_inf);
        global::argv(6).read(uMax);
        global::argv(7).read(outDir);
    } catch(PlbIOException& exception) {
        pcout << exception.what() << endl;
        pcout << "Command line arguments:\n";
        pcout << "1 : d_part\n";
        pcout << "2 : N per particle diameter\n";
        pcout << "3 : particle volume fraction\n";
        pcout << "4 : nu_fluid\n";
        pcout << "5 : estimated v_inf\n";
        pcout << "6 : uMax\n";
        pcout << "7 : outDir\n";
        exit(1);
    }

    std::string lbOutDir(outDir), demOutDir(outDir);
    lbOutDir.append("tmp/");
    demOutDir.append("post/");
    global::directories().setOutputDir(lbOutDir);

    const T rho_f = 1000;

    LiggghtsCouplingWrapper wrapper(argv,global::mpi().getGlobalCommunicator());

    // particle size and volume fraction are handed over to LIGGGHTS
    // as variables (see LIGGGHTS docu for details)
    wrapper.setVariable("r_part",d_part/2);
    wrapper.setVariable("v_frac",v_frac);

    wrapper.execFile("in.lbdem");


    T g = 9.81;

    const T lx = 1., ly = 1., lz = 2.;


    T r_ = d_part/2.;
    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;


//.........这里部分代码省略.........