毕业了——课题代码开源（二）使用Palabos的单相流多孔渗流代码

小咪咪 2023-02-16 09:27 31阅读 0赞

# 简介 #

这是我依据[Palabos案例解析（一）permeability.cpp案例][Palabos_permeability.cpp]的原始案例修改的适合我自己的代码，其实没咋动，主要就是改一下最后速度的统计的标准和数据输出的一些小东西，其实和源代码几乎没啥区别。

# 代码 #

#include "palabos3D.h"
    #include "palabos3D.hh"
    
    #include <vector>
    #include <cmath>
    #include <cstdlib>
    
    using namespace plb;
    
    typedef double T;
    #define DESCRIPTOR descriptors::D3Q19Descriptor
    
    // This function object returns a zero velocity, and a pressure which decreases
    //   linearly in z-direction. It is used to initialize the particle populations.
    class PressureGradient {
        
    public:
        PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
        {
         }
        void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
        {
        
            velocity.resetToZero();
            density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;
        }
    private:
        T deltaP;
        plint nx;
    };
    
    void readGeometry(std::string fNameIn, std::string fNameOut, MultiScalarField3D<int>& geometry)
    {
        
        const plint nx = geometry.getNx();
        const plint ny = geometry.getNy();
        const plint nz = geometry.getNz();
    
        Box3D sliceBox(0,0, 0,ny-1, 0,nz-1);
        std::auto_ptr<MultiScalarField3D<int> > slice = generateMultiScalarField<int>(geometry, sliceBox);
        plb_ifstream geometryFile(fNameIn.c_str());
        for (plint iX=0; iX<nx; ++iX) {
        
            if (!geometryFile.is_open()) {
        
                pcout << "Error: could not open geometry file " << fNameIn << std::endl;
                exit(EXIT_FAILURE);
            }
            geometryFile >> *slice;
            copy(*slice, slice->getBoundingBox(), geometry, Box3D(iX,iX, 0,ny-1, 0,nz-1));
        }
    
        {
        
            VtkImageOutput3D<T> vtkOut("porousMedium", 1.0);
            vtkOut.writeData<float>(*copyConvert<int,T>(geometry, geometry.getBoundingBox()), "tag", 1.0);
        }
    
        {
        
            /*std::auto_ptr<MultiScalarField3D<T> > floatTags = copyConvert<int,T>(geometry, geometry.getBoundingBox());
            std::vector<T> isoLevels;
            isoLevels.push_back(0.5);
            typedef TriangleSet<T>::Triangle Triangle;
            std::vector<Triangle> triangles;
            Box3D domain = floatTags->getBoundingBox().enlarge(-1);
            domain.x0++;
            domain.x1--;
            isoSurfaceMarchingCube(triangles, *floatTags, isoLevels, domain);
            TriangleSet<T> set(triangles);
            std::string outDir = fNameOut + "/";
            set.writeBinarySTL(outDir + "porousMedium.stl");*/
        }
    }
    
    void porousMediaSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
            OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition,
            MultiScalarField3D<int>& geometry, T deltaP)
    {
        
        const plint nx = lattice.getNx();
        const plint ny = lattice.getNy();
        const plint nz = lattice.getNz();
    
        pcout << "Definition of inlet/outlet." << std::endl;
        Box3D inlet (1,nx-2, 1,ny-2, 0,0);
        boundaryCondition->addPressureBoundary0N(inlet, lattice);
        setBoundaryDensity(lattice, inlet, (T) 1.);
    
        Box3D outlet(1,nx-2, 1,ny-2, nz-1,nz-1);
        boundaryCondition->addPressureBoundary0P(outlet, lattice);
        setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2);
    
        pcout << "Definition of the geometry." << std::endl;
        // Where "geometry" evaluates to 1, use bounce-back.
        defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
        // Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
        defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);
    
        Box3D X0Plane(0,0,0,ny-1,0,nz-1);
        Box3D XnPlane(nx-1,nx-1,0,ny-1,0,nz-1);
        Box3D Y0Plane(0,nx-1,0,0,0,nz-1);
        Box3D YnPlane(0,nx-1,ny-1,ny-1,0,nz-1);
        defineDynamics(lattice, X0Plane, new BounceBack<T,DESCRIPTOR>());
        defineDynamics(lattice, XnPlane, new BounceBack<T,DESCRIPTOR>());
        defineDynamics(lattice, Y0Plane, new BounceBack<T,DESCRIPTOR>());
        defineDynamics(lattice, YnPlane, new BounceBack<T,DESCRIPTOR>());
    
        pcout << "Initilization of rho and u." << std::endl;
        initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) );
    
        lattice.initialize();
        delete boundaryCondition;
    }
    
    void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
    {
        
        const plint nx = lattice.getNx();
        const plint ny = lattice.getNy();
        const plint nz = lattice.getNz();
    
        const plint imSize = 600;
        ImageWriter<T> imageWriter("leeloo");
    
        // Write velocity-norm at x=0.
        imageWriter.writeScaledGif(createFileName("ux_inlet", iter, 6),
                *computeVelocityNorm(lattice, Box3D(0,0, 0,ny-1, 0,nz-1)),
                imSize, imSize );
    
        // Write velocity-norm at x=nx/2.
        imageWriter.writeScaledGif(createFileName("ux_half", iter, 6),
                *computeVelocityNorm(lattice, Box3D(nx/2,nx/2, 0,ny-1, 0,nz-1)),
                imSize, imSize );
    }
    
    void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
    {
        
        VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), 1.);
        vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", 1.);
        vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", 1.);
    }
    
    T computePermeability(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box3D domain )
    {
        
        pcout << "Computing the permeability." << std::endl;
    
        // Compute only the z-direction of the velocity (direction of the flow).
        plint zComponent = 2;
        plint nz = lattice.getNz();
    
        T meanU = computeAverage(*computeVelocityComponent(lattice, domain, zComponent));
    
        pcout << "Average velocity     = " << meanU                         << std::endl;
        pcout << "Lattice viscosity nu = " << nu                            << std::endl;
        pcout << "Grad P               = " << deltaP/(T)(nz-1)              << std::endl;
        pcout << "Permeability         = " << nu*meanU / (deltaP/(T)(nz-1)) << std::endl;
    
        return meanU;
    }
    
    int main(int argc, char **argv)
    {
        
        plbInit(&argc, &argv);
    
        if (argc!=7) {
        
            pcout << "Error missing some input parameter\n";
            pcout << "The structure is :\n";
            pcout << "1. Input file name.\n";
            pcout << "2. Output directory name.\n";
            pcout << "3. number of cells in X direction.\n";
            pcout << "4. number of cells in Y direction.\n";
            pcout << "5. number of cells in Z direction.\n";
            pcout << "6. Delta P .\n";
            pcout << "Example: " << argv[0] << " palabos.dat tmp/ 201 201 201 0.00005\n";
            exit (EXIT_FAILURE);
        }
        std::string fNameIn  = argv[1];
        std::string fNameOut = argv[2];
    
        const plint nx = atoi(argv[3]);
        const plint ny = atoi(argv[4]);
        const plint nz = atoi(argv[5]);
        const T deltaP = atof(argv[6]);
    
        global::directories().setOutputDir(fNameOut+"/");
    
        const T omega = 1.0;
        const T nu    = ((T)1/omega- (T)0.5)/DESCRIPTOR<T>::invCs2;
    
        pcout << "Creation of the lattice." << std::endl;
        MultiBlockLattice3D<T,DESCRIPTOR> lattice(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega));
        // Switch off periodicity.
        lattice.periodicity().toggleAll(false);
    
        pcout << "Reading the geometry file." << std::endl;
        MultiScalarField3D<int> geometry(nx,ny,nz);
        readGeometry(fNameIn, fNameOut, geometry);
    
        pcout << "nu = " << nu << std::endl;
        pcout << "deltaP = " << deltaP << std::endl;
        pcout << "omega = " << omega << std::endl;
        pcout << "nx = " << lattice.getNx() << std::endl;
        pcout << "ny = " << lattice.getNy() << std::endl;
        pcout << "nz = " << lattice.getNz() << std::endl;
    
        porousMediaSetup(lattice, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);
    
        // The value-tracer is used to stop the simulation once is has converged.
        // 1st parameter:velocity
        // 2nd parameter:size
        // 3rd parameters:threshold
        // 1st and second parameters ae used for the length of the time average (size/velocity)
        util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);
    
        pcout << "Simulation begins" << std::endl;
        plint iT=0;
    
        const plint maxT = 30000;
        for (;iT<maxT; ++iT) {
        
            if (iT % 20 == 0) {
        
                pcout << "Iteration " << iT << std::endl;
            }
            /*if (iT % 500 == 0 && iT>0) {
                writeGifs(lattice,iT);
            }*/
    
            lattice.collideAndStream();
            converge.takeValue(getStoredAverageEnergy(lattice),true);
    
            if (converge.hasConverged()) {
        
                break;
            }
        }
    
        pcout << "End of simulation at iteration " << iT << std::endl;
    
        pcout << "Permeability:" << std::endl << std::endl;
        computePermeability(lattice, nu, deltaP, lattice.getBoundingBox());
        pcout << std::endl;
    
        pcout << "Writing VTK file ..." << std::endl << std::endl;
        writeVTK(lattice,iT);
        pcout << "Finished!" << std::endl << std::endl;
    
        return 0;
    }

# 说明 #

其实这个是没啥可说的，对于参数，我并没有特别的调过，还需要你自己依照自己课题的对应尺寸去调，图片就不放了，跟那个渗透的案例差不多，想看图的可以看那个的图。

# 附加题 #

## 代码 ##

#include "palabos2D.h"
    #include "palabos2D.hh"
    
    #include <vector>
    #include <cmath>
    #include <cstdlib>
    
    using namespace plb;
    
    typedef double T;
    #define DESCRIPTOR descriptors::D2Q9Descriptor
    
    // This function object returns a zero velocity, and a pressure which decreases
    //   linearly in z-direction. It is used to initialize the particle populations.
    class PressureGradient {
        
    public:
        PressureGradient(T deltaP_, plint ny_) : deltaP(deltaP_), ny(ny_)
        {
         }
        void operator() (plint iX, plint iY, T& density, Array<T,2>& velocity) const
        {
        
            velocity.resetToZero();
            density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(ny-1) * (T)(ny-1-iY);
        }
    private:
        T deltaP;
        plint ny;
    };
    
    void readGeometry(std::string fNameIn, std::string fNameOut, MultiScalarField2D<int>& geometry)
    {
        
        plb_ifstream geometryFile(fNameIn.c_str());
        if (!geometryFile.is_open()) {
        
            pcout << "Error: could not open geometry file " << fNameIn << std::endl;
            exit(EXIT_FAILURE);
        }
        geometryFile >> geometry;
    
        {
        
            VtkImageOutput2D<T> vtkOut("porousMedium", 1.0);
            vtkOut.writeData<float>(*copyConvert<int,T>(geometry, geometry.getBoundingBox()), "tag", 1.0);
        }
    }
    
    void porousMediaSetup(MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
            OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition,
            MultiScalarField2D<int>& geometry, T deltaP)
    {
        
        const plint nx = lattice.getNx();
        const plint ny = lattice.getNy();
    
        pcout << "Definition of inlet/outlet." << std::endl;
        Box2D inlet (1,nx-2, ny-1,ny-1);
        boundaryCondition->addPressureBoundary1P(inlet, lattice);
        setBoundaryDensity(lattice, inlet, (T) 1.);
    
        Box2D outlet(1,nx-2, 0,0);
        boundaryCondition->addPressureBoundary1N(outlet, lattice);
        setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2);
    
        pcout << "Definition of the geometry." << std::endl;
        // Where "geometry" evaluates to 1, use bounce-back.
        defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
        // Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
        defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);
    
        pcout << "Initilization of rho and u." << std::endl;
        initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, ny) );
    
        lattice.initialize();
        delete boundaryCondition;
    }
    
    void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
    {
        
        const plint nx = lattice.getNx();
        const plint ny = lattice.getNy();
    
        const plint imSize = 600;
        ImageWriter<T> imageWriter("leeloo");
    
        imageWriter.writeScaledGif(createFileName("u", iter, 6),
                *computeVelocityNorm(lattice, Box2D(0,nx-1, 0,ny-1)),
                imSize, imSize );
    }
    
    void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
    {
        
        VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 6), 1.);
        vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", 1.);
        vtkOut.writeData<2,float>(*computeVelocity(lattice), "velocity", 1.);
    }
    
    T computePermeability(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box2D domain )
    {
        
        pcout << "Computing the permeability." << std::endl;
    
        // Compute only the y-direction of the velocity (direction of the flow).
        plint yComponent = 1;
        plint ny = lattice.getNy();
    
        T meanU = -computeAverage(*computeVelocityComponent(lattice, domain, yComponent));
    
        pcout << "Average velocity     = " << meanU                         << std::endl;
        pcout << "Lattice viscosity nu = " << nu                            << std::endl;
        pcout << "Grad P               = " << deltaP/(T)(ny-1)              << std::endl;
        pcout << "Permeability         = " << nu*meanU / (deltaP/(T)(ny-1)) << std::endl;
    
        return meanU;
    }
    
    int main(int argc, char **argv)
    {
        
        plbInit(&argc, &argv);
    
        if (argc!=6) {
        
            pcout << "Error missing some input parameter\n";
            pcout << "The structure is :\n";
            pcout << "1. Input file name.\n";
            pcout << "2. Output directory name.\n";
            pcout << "3. number of cells in X direction.\n";
            pcout << "4. number of cells in Y direction.\n";
            pcout << "5. Delta P .\n";
            pcout << "Example: " << argv[0] << " palabos.dat tmp/ 201 201 0.00005\n";
            exit (EXIT_FAILURE);
        }
        std::string fNameIn  = argv[1];
        std::string fNameOut = argv[2];
    
        const plint nx = atoi(argv[3]);
        const plint ny = atoi(argv[4]);
        const T deltaP = atof(argv[5]);
    
        global::directories().setOutputDir(fNameOut+"/");
    
        const T omega = 1.0;
        const T nu    = ((T)1/omega- (T)0.5)/DESCRIPTOR<T>::invCs2;
    
        pcout << "Creation of the lattice." << std::endl;
        MultiBlockLattice2D<T,DESCRIPTOR> lattice(nx,ny, new BGKdynamics<T,DESCRIPTOR>(omega));
        // Switch off periodicity.
        lattice.periodicity().toggleAll(false);
    
        pcout << "Reading the geometry file." << std::endl;
        MultiScalarField2D<int> geometry(nx,ny);
        readGeometry(fNameIn, fNameOut, geometry);
    
        pcout << "nu = " << nu << std::endl;
        pcout << "deltaP = " << deltaP << std::endl;
        pcout << "omega = " << omega << std::endl;
        pcout << "nx = " << lattice.getNx() << std::endl;
        pcout << "ny = " << lattice.getNy() << std::endl;
    
        porousMediaSetup(lattice, createLocalBoundaryCondition2D<T,DESCRIPTOR>(), geometry, deltaP);
    
        // The value-tracer is used to stop the simulation once is has converged.
        // 1st parameter:velocity
        // 2nd parameter:size
        // 3rd parameters:threshold
        // 1st and second parameters ae used for the length of the time average (size/velocity)
        util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);
    
        pcout << "Simulation begins" << std::endl;
        plint iT=0;
    
        const plint maxT = 100000;
        for (;iT<maxT; ++iT) {
        
            if (iT % 20 == 0) {
        
                pcout << "Iteration " << iT << std::endl;
            }
            if (iT % 50 == 0 && iT>0) {
        
                writeGifs(lattice,iT);
            }
    
            lattice.collideAndStream();
            converge.takeValue(getStoredAverageEnergy(lattice),true);
    
            if (converge.hasConverged()) {
        
                break;
            }
        }
    
        pcout << "End of simulation at iteration " << iT << std::endl;
    
        pcout << "Permeability:" << std::endl << std::endl;
        computePermeability(lattice, nu, deltaP, lattice.getBoundingBox());
        pcout << std::endl;
    
        pcout << "Writing VTK file ..." << std::endl << std::endl;
        writeVTK(lattice,iT);
        pcout << "Finished!" << std::endl << std::endl;
    
        return 0;
    }

## 说明 ##

这是二维的渗透，仅仅是拿三维直接改的版本，对于参数，我也是并没有特别的细究，对于此程序，可以和上面三维的和[原始渗透案例][Palabos_permeability.cpp]进行一下对比，你就明白在palabos中如何三维转二维和二维转三维。

# 完 #

[Palabos_permeability.cpp]: https://blog.csdn.net/qq_28632981/article/details/103756787