LBM：PalaBos-OpenLB 提前体验自由表面流

xhbj

From我的博客:http://blog.sina.com.cn/vonkarman
LBM PalaBos 新版本将在1周内发布, 敬请期待！
以下是新版本中的自由表面流的在例子（流体为水）。 敬请享用……

计算图像见： http://blog.sina.com.cn/vonkarman

#include "palabos3D.h"
#include "palabos3D.hh"

using namespace plb;
using namespace std;

#define DESCRIPTOR descriptors::ForcedD3Q19Descriptor
typedef double T;


// Smagorinsky constant for LES model.
const T cSmago = 0.14;

// Physical dimensions of the system (in meters).
const T lx = 3.22;
const T ly = 1.0;
const T lz = 1.0;

const T iniRhoFluid = T(1);
const T rhoEmpty    = T(1);
Array<T,3> forceOrientation(T(),T(),(T)1);
   
plint writeImagesIter   = 100;
plint getStatisticsIter = 20;

plint maxIter;
plint N;
plint nx, ny, nz;
T delta_t, delta_x;
Array<T,3> externalForce;
T nuPhys, nuLB, tau, omega;

std::string outDir;
plint obstacleCenterXYplane, obstacleLength, obstacleWidth, obstacleHeight, beginWaterReservoir, waterReservoirHeight;
plint waterLevelOne, waterLevelTwo, waterLevelThree, waterLevelFour;

void setupParameters() {
    delta_x = T(1)/N;
    nx = (plint) (lx*N);
    ny = (plint) (ly*N);
    nz = (plint) (lz*N);

    // Gravity in lattice units.
    externalForce = Array<T,3>(0., 0., (-9.8* delta_t * delta_t)/delta_x);
    tau            = (nuPhys*DESCRIPTOR<T>::invCs2*delta_t)/(delta_x*delta_x) + 0.5;
    omega          = 1./tau;   
    nuLB           = (tau-0.5)*DESCRIPTOR<T>::cs2; // Viscosity in lattice units.
   
    obstacleCenterXYplane = util::roundToInt(0.744*N);
    obstacleLength        = util::roundToInt(0.403*N);
    obstacleWidth         = util::roundToInt(0.161*N);
    obstacleHeight        = util::roundToInt(0.161*N);
    beginWaterReservoir   = util::roundToInt((0.744+1.248)*N);
    waterReservoirHeight  = util::roundToInt(0.55*N);
   
    waterLevelOne   = util::roundToInt(0.496*N);
    waterLevelTwo   = util::roundToInt(2.*0.496*N);
    waterLevelThree = util::roundToInt(3.*0.496*N);
    waterLevelFour  = util::roundToInt((3.*0.496 + 1.150)*N);
}

// Specifies the initial condition for the fluid (each cell is assigned the
// flag "fluid", "empty", or "wall").
int initialFluidFlags(plint iX, plint iY, plint iZ) {
    // Place an obstacle on the left end, which is hit by the fluid.
    bool insideObstacle =
        iX >= obstacleCenterXYplane-obstacleWidth/2 &&
        iX <= obstacleCenterXYplane+obstacleWidth/2 &&
        iY >= ny/2-obstacleLength/2 &&
        iY <= ny/2+obstacleLength/2 &&
        iZ <= obstacleHeight+1;
   
    if (insideObstacle) {
        return freeSurfaceFlag::wall;
    }
    else if (iX >= beginWaterReservoir && iZ <= waterReservoirHeight) {
        return freeSurfaceFlag::fluid;
    }
    else {
        return freeSurfaceFlag::empty;
    }
}

void writeResults(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, MultiScalarField3D<T>& volumeFraction, plint iT)
{
    static const plint nx = lattice.getNx();
    static const plint ny = lattice.getNy();
    static const plint nz = lattice.getNz();
    Box3D slice(0, nx-1, ny/2, ny/2, 0, nz-1);
    ImageWriter<T> imageWriter("leeloo");
    imageWriter.writeScaledPpm(createFileName("u", iT, 6),
                               *computeVelocityNorm(lattice, slice));

    imageWriter.writeScaledPpm(createFileName("rho", iT, 6),
                               *computeDensity(lattice, slice));
                  
    imageWriter.writeScaledPpm(createFileName("volumeFraction", iT, 6), *extractSubDomain(volumeFraction, slice));

    // Use a marching-cube algorithm to reconstruct the free surface and write an STL file.
    std::vector<T> isoLevels;
    isoLevels.push_back((T) 0.5);
    typedef TriangleSet<T>::Triangle Triangle;
    std::vector<Triangle> triangles;
    isoSurfaceMarchingCube(triangles, volumeFraction, isoLevels, volumeFraction.getBoundingBox().enlarge(-5));
    TriangleSet<T>(triangles).writeAsciiSTL(createFileName(outDir+"/interface", iT, 6)+".stl");
}

void writeStatistics(FreeSurfaceFields3D<T,DESCRIPTOR>& fields) {
    pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << endl;
    T averageMass = freeSurfaceAverageMass<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());
    pcout << "Average Mass: " << averageMass  << endl;
    T averageDensity =freeSurfaceAverageDensity<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());
    pcout << "Average Density: " << setprecision(12) << averageDensity  << endl;

    T averageVolumeFraction = freeSurfaceAverageVolumeFraction<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());
    pcout << "Average Volume-Fraction: " << setprecision(12) << averageVolumeFraction  << endl;

    T heightAtOne = getAverageHeightAtXY<T,DESCRIPTOR> (
            fields.freeSurfaceArgs, N,
            Box3D(waterLevelOne, waterLevelOne, ny/2, ny/2, 1, nz-1) );
    pcout << "Height at pos 1: " << setprecision(12) << heightAtOne  << endl;

    T heightAtTwo = getAverageHeightAtXY<T,DESCRIPTOR> (
            fields.freeSurfaceArgs, N,
            Box3D(waterLevelTwo,waterLevelTwo,ny/2,ny/2,1,nz-1) );
    pcout << "Height at pos 2: " << setprecision(12) << heightAtTwo  << endl;
   
    T heightAtThree = getAverageHeightAtXY<T,DESCRIPTOR> (
            fields.freeSurfaceArgs, N,
            Box3D(waterLevelThree,waterLevelThree,ny/2,ny/2,1,nz-1) );
    pcout << "Height at pos 3: " << setprecision(12) << heightAtThree  << endl;
   
    T heightAtFour = getAverageHeightAtXY<T,DESCRIPTOR> (
            fields.freeSurfaceArgs, N,
            Box3D(waterLevelFour,waterLevelFour, ny/2,ny/2, 1,nz-1) );
    pcout << "Height at pos 4: " << setprecision(12) << heightAtFour  << endl;
    pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << endl;
}


int main(int argc, char **argv)
{
    plbInit(&argc, &argv);
    global::directories().setInputDir("./");
        
    if (global::argc() != 6) {
        pcout << "Error missing some input parameter\n";
    }

    try {
        global::argv(1).read(outDir);
        global::directories().setOutputDir(outDir+"/");

        global::argv(2).read(nuPhys);
        global::argv(3).read(N);
        global::argv(4).read(delta_t);
        global::argv(5).read(maxIter);
    }
    catch(PlbIOException& except) {
        pcout << except.what() << std::endl;
        pcout << "The parameters for this program are :\n";
        pcout << "1. Output directory name.\n";
        pcout << "2. kinematic viscosity in physical Units (m^2/s) .\n";
        pcout << "3. number of lattice nodes for lz .\n";
        pcout << "4. delta_t .\n";
        pcout << "5. maxIter .\n";
        pcout << "Reasonable parameters on a desktop computer are: " << (std::string)global::argv(0) << " tmp 1.e-5 40 1.e-3 80000\n";
        pcout << "Reasonable parameters on a parallel machine are: " << (std::string)global::argv(0) << " tmp 1.e-6 100 1.e-4 80000\n";
        exit (EXIT_FAILURE);
    }
   
    setupParameters();
   
    pcout << "delta_t= " << delta_t << endl;
    pcout << "delta_x= " << delta_x << endl;
    pcout << "delta_t*delta_t/delta_x= " << delta_t*delta_t/delta_x << endl;
    pcout << "externalForce= " << externalForce[2] << endl;
    pcout << "relaxation time= " << tau << endl;
    pcout << "omega= " << omega << endl;
    pcout << "kinematic viscosity physical units = " << nuPhys << endl;
    pcout << "kinematic viscosity lattice units= " << nuLB << endl;
   
    global::timer("initialization").start();
   

    SparseBlockStructure3D blockStructure(createRegularDistribution3D(nx, ny, nz));

    Dynamics<T,DESCRIPTOR>* dynamics
        = new SmagorinskyBGKdynamics<T,DESCRIPTOR>(omega, cSmago);

    FreeSurfaceFields3D<T,DESCRIPTOR> fields(blockStructure, dynamics->clone(), iniRhoFluid, rhoEmpty, externalForce);

    // Set all outer-wall cells to "wall" (here, bulk-cells are also set to "wall", but it
    // doesn't matter, because they are overwritten on the next line).
    setToConstant(fields.flag, fields.flag.getBoundingBox(), (int)freeSurfaceFlag::wall);
    // In the bulk (all except outer wall layer), initialize the flags as specified by
    // the function "initialFluidFlags".
    setToFunction(fields.flag, fields.flag.getBoundingBox().enlarge(-1), initialFluidFlags);
   
    fields.defaultInitialize();

    pcout << "Time spent for setting up lattices: "
          << global::timer("initialization").stop() << endl;
    T lastIterationTime = T();

    for (plint iT = 0; iT <= maxIter; ++iT) {
        global::timer("iteration").restart();
        
        if (iT % getStatisticsIter==0) {
            pcout << endl;
            pcout << "ITERATION = " << iT << endl;
            pcout << "Time of last iteration is " << lastIterationTime << " seconds" << endl;
            writeStatistics(fields);
        }
        
        // This includes the collision-streaming cycle, plus all free-surface operations.
        fields.lattice.executeInternalProcessors();

        if (iT % writeImagesIter == 0) {
            global::timer("images").start();
            writeResults(fields.lattice, fields.volumeFraction, iT);
            pcout << "Total time spent for writing images: "
                << global::timer("images").stop() << endl;
        }                           
        lastIterationTime = global::timer("iteration").stop();
    }
}