Skip to content

Pressure Boundary Condition

  • This topic has 2 replies, 1 voice, and was last updated 8 years ago by ivan.
Viewing 3 posts - 1 through 3 (of 3 total)
  • Author
    Posts
  • #1863
    ivan
    Member

    Okay, so here is one of the combination I’ve been trying recently with which I simply get a lot of “nan” results. The bold lines are the ones related to the pressure boundary condition I am trying to include and without them the simulation works.

    Code:
    #include “olb2D.h”
    #include “olb2D.hh”
    #include <cstdlib>
    #include <iostream>

    using namespace olb;
    using namespace olb::descriptors;
    using namespace olb::graphics;
    using namespace std;

    typedef double T;
    #define DESCRIPTOR ShanChenDynOmegaForcedD2Q9Descriptor

    // Parameters for the simulation setup
    const int nx = 100;
    const int ny = 250;
    const T centerDropletX = (nx-1)/2;
    const T centerDropletY = 20.;
    const T radiusDroplet = 10.;
    const T wallThickness =4.;
    const int maxIter = 6900;

    /// Stores geometry information in form of material numbers
    void prepareGeometry(SuperGeometry2D<T>& superGeometry) {

    OstreamManager clout(std::cout,”prepareGeometry”);
    clout << “Prepare Geometry …” << std::endl;

    // Sets material number for fluid and boundary
    superGeometry.rename(0,1);

    std::vector<T> origin1(2,T());
    origin1[0] = 0;
    origin1[1] = 0;

    std::vector<T> origin2(2,T());
    origin2[0] = 0;
    origin2[1] = ny-1;

    std::vector<T> origin3(2,T());
    origin3[0] = 0;
    origin3[1] = wallThickness;

    std::vector<T> origin4(2,T());
    origin4[0] = nx-wallThickness-1;
    origin4[1] = wallThickness;

    std::vector<T> extend1(2,T());
    extend1[0] = nx;
    extend1[1] = wallThickness;

    std::vector<T> extend2(2,T());
    extend2[0] = wallThickness;
    extend2[1] = ny-wallThickness;

    std::vector<T>extend3(2,T());
    extend3[0] = nx;
    extend3[1] = 1;

    Vector<T,2> centerDroplet(centerDropletX,centerDropletY);

    IndicatorCuboid2D<T> bottom(extend1, origin1);
    IndicatorCuboid2D<T> top(extend3, origin2);
    IndicatorCuboid2D<T> leftwall(extend2, origin3);
    IndicatorCuboid2D<T> rightwall(extend2, origin4);
    IndicatorCircle2D<T> droplet(centerDroplet, radiusDroplet);

    superGeometry.rename(1,2,droplet);
    superGeometry.rename(1,3,bottom);
    superGeometry.rename(1,3,leftwall);
    superGeometry.rename(1,3,rightwall);
    superGeometry.rename(1,4,top);

    /// Removes all not needed boundary voxels inside the surface
    superGeometry.innerClean();
    superGeometry.checkForErrors();

    superGeometry.print();

    clout << “Prepare Geometry … OK” << std::endl;
    }

    void prepareLattice(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    Dynamics<T, DESCRIPTOR>& bulkDynamics1,
    Dynamics<T, DESCRIPTOR>& bulkDynamics2,
    Dynamics<T, DESCRIPTOR>& bounceBackRho0,
    Dynamics<T, DESCRIPTOR>& bounceBackRho1,
    [b] sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeOne,
    sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeOne,
    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeTwo,
    sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeTwo,[/b]
    SuperGeometry2D<T>& superGeometry)
    {

    OstreamManager clout(std::cout,”prepareLattice”);
    clout << “Prepare Lattice …” << std::endl;

    [b] const T omega1 = 1.0;
    const T omega2 = 1.0;[/b]

    /// define lattice Dynamics
    sLatticeOne.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
    sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
    sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
    sLatticeOne.defineDynamics(superGeometry, 3, &bounceBackRho1);
    sLatticeOne.defineDynamics(superGeometry, 4, &bulkDynamics1);

    [b] sBoundaryConditionLatticeOne.addPressureBoundary(superGeometry, 4, omega1);[/b]

    sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
    sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
    sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
    sLatticeTwo.defineDynamics(superGeometry, 3, &bounceBackRho1);
    sLatticeTwo.defineDynamics(superGeometry, 4, &bulkDynamics2);
    [b]
    sBoundaryConditionLatticeTwo.addPressureBoundary(superGeometry, 4, omega2);[/b]

    clout << “Prepare Lattice … OK” << std::endl;
    }

    void setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    T force, int iT, SuperGeometry2D<T>& superGeometry) {

    if(iT==0) {

    AnalyticalConst2D<T,T> noise(4.e-2);
    std::vector<T> v(2,T());
    AnalyticalConst2D<T,T> zeroV(v);
    AnalyticalConst2D<T,T> zero(0.);
    AnalyticalLinear2D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.98+force*ny*DESCRIPTOR<T>::invCs2);
    AnalyticalConst2D<T,T> onePlus(0.98+force*ny/2.*DESCRIPTOR<T>::invCs2);
    AnalyticalRandom2D<T,T> random;
    AnalyticalIdentity2D<T,T> randomOne(random*noise+one);
    AnalyticalIdentity2D<T,T> randomPlus(random*noise+onePlus);
    std::vector<T> F(2,T());
    F[1] = -force;
    AnalyticalConst2D<T,T> f(F);

    /// for each material set the defineRhou and the Equilibrium

    sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 1,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 2,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeOne.defineRhoU(superGeometry, 4, zero, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 4, zero, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 4,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);

    sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);

    sLatticeTwo.defineRhoU(superGeometry, 4, randomPlus, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 4, randomPlus, zeroV);

    /// Make the lattice ready for simulation
    sLatticeOne.initialize();
    sLatticeTwo.initialize();
    }
    }

    void getResults(SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeOne, int iT,
    SuperGeometry2D<T>& superGeometry, Timer<T>& timer) {

    OstreamManager clout(std::cout,”getResults”);
    SuperVTKwriter2D<T> vtkWriter(“rayleighTaylor2dsLatticeOne”);

    const int vtkIter = 100;
    const int statIter = 10;

    if (iT==0) {
    /// Writes the geometry, cuboid no. and rank no. as vti file for visualization
    SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLatticeOne, superGeometry);
    SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLatticeOne);
    SuperLatticeRank2D<T, DESCRIPTOR> rank(sLatticeOne);
    vtkWriter.write(geometry);
    vtkWriter.write(cuboid);
    vtkWriter.write(rank);
    vtkWriter.createMasterFile();
    }

    // Get statistics
    if (iT%statIter==0 && iT > 0) {
    // Timer console output
    timer.update( iT );
    timer.printStep();

    clout << “averageRhoFluidOne=” << sLatticeOne.getStatistics().getAverageRho();
    clout << “; averageRhoFluidTwo=” << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
    }

    /// Writes the VTK files
    if (iT%vtkIter==0) {
    clout << “Writing VTK …” << std::endl;
    SuperLatticeVelocity2D<T, DESCRIPTOR> velocity(sLatticeOne);
    SuperLatticeDensity2D<T, DESCRIPTOR> density(sLatticeOne);
    vtkWriter.addFunctor( velocity );
    vtkWriter.addFunctor( density );
    vtkWriter.write(iT);

    BlockLatticeReduction2D<T, DESCRIPTOR> planeReduction(density);
    BlockGifWriter<T> gifWriter;
    gifWriter.write(planeReduction, iT, “density”);

    clout << “Writing VTK … OK” << std::endl;
    }
    }

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

    /// === 1st Step: Initialization ===

    olbInit(&argc, &argv);
    singleton::directories().setOutputDir(“./tmp/”);
    OstreamManager clout(std::cout,”main”);

    const T omega1 = 1.0;
    const T omega2 = 1.0;
    const T G = 3.;
    T force = -30./(T)ny/(T)ny;

    /// === 2nd Step: Prepare Geometry ===
    /// Instantiation of a cuboidGeometry with weights

    #ifdef PARALLEL_MODE_MPI
    CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, singleton::mpi().getSize());
    #else
    CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, 1);
    #endif

    cGeometry.setPeriodicity(false, false);

    HeuristicLoadBalancer<T> loadBalancer(cGeometry);

    SuperGeometry2D<T> superGeometry(cGeometry, loadBalancer, 2);

    prepareGeometry(superGeometry);

    /// === 3rd Step: Prepare Lattice ===

    SuperLattice2D<T, DESCRIPTOR> sLatticeOne(superGeometry);
    SuperLattice2D<T, DESCRIPTOR> sLatticeTwo(superGeometry);

    ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
    omega1, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );
    ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics2 (
    omega2, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );

    // A bounce-back node with fictitious density 1,
    // which is experienced by the partner fluid
    BounceBack<T, DESCRIPTOR> bounceBackRho1( 1. );
    // A bounce-back node with fictitious density 0,
    // which is experienced by the partner fluid
    BounceBack<T, DESCRIPTOR> bounceBackRho0( 0. );

    std::vector<T> rho0;
    rho0.push_back(1);
    rho0.push_back(1);
    PsiEqualsRho<T,T> interactionPotential;
    ShanChenForcedGenerator2D<T,DESCRIPTOR> coupling(G,rho0,interactionPotential);

    [b] // Pressure boundary
    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeOne(sLatticeOne);
    createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeOne);
    sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeOne(sLatticeOne);

    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeTwo(sLatticeTwo);
    createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeTwo);
    sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeTwo(sLatticeTwo);[/b]

    sLatticeOne.addLatticeCoupling(superGeometry, 1, coupling, sLatticeTwo);
    sLatticeOne.addLatticeCoupling(superGeometry, 2, coupling, sLatticeTwo);

    prepareLattice(sLatticeOne, sLatticeTwo, bulkDynamics1, bulkDynamics2,
    bounceBackRho0, bounceBackRho1,[b] sBoundaryConditionLatticeOne, sOffBoundaryConditionLatticeOne, sBoundaryConditionLatticeTwo, sOffBoundaryConditionLatticeTwo,[/b] superGeometry);

    /// === 4th Step: Main Loop with Timer ===
    int iT = 0;
    clout << “starting simulation…” << endl;
    Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
    timer.start();

    for (iT=0; iT<maxIter; ++iT) {

    /// === 5th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues(sLatticeOne, sLatticeTwo, force, iT, superGeometry);

    /// === 6th Step: Collide and Stream Execution ===
    sLatticeOne.collideAndStream();
    sLatticeTwo.collideAndStream();

    sLatticeOne.communicate();
    sLatticeTwo.communicate();

    sLatticeOne.executeCoupling();

    /// === 7th Step: Computation and Output of the Results ===
    getResults(sLatticeTwo, sLatticeOne, iT, superGeometry, timer);
    }

    timer.stop();
    timer.printSummary();
    }

    #2453
    ivan
    Member

    Okay, so here is one of the combination I’ve been trying recently with which I simply get a lot of “nan” results. The bold lines are the ones related to the pressure boundary condition I am trying to include and without them the simulation works.

    Code:
    #include “olb2D.h”
    #include “olb2D.hh”
    #include <cstdlib>
    #include <iostream>

    using namespace olb;
    using namespace olb::descriptors;
    using namespace olb::graphics;
    using namespace std;

    typedef double T;
    #define DESCRIPTOR ShanChenDynOmegaForcedD2Q9Descriptor

    // Parameters for the simulation setup
    const int nx = 100;
    const int ny = 250;
    const T centerDropletX = (nx-1)/2;
    const T centerDropletY = 20.;
    const T radiusDroplet = 10.;
    const T wallThickness =4.;
    const int maxIter = 6900;

    /// Stores geometry information in form of material numbers
    void prepareGeometry(SuperGeometry2D<T>& superGeometry) {

    OstreamManager clout(std::cout,”prepareGeometry”);
    clout << “Prepare Geometry …” << std::endl;

    // Sets material number for fluid and boundary
    superGeometry.rename(0,1);

    std::vector<T> origin1(2,T());
    origin1[0] = 0;
    origin1[1] = 0;

    std::vector<T> origin2(2,T());
    origin2[0] = 0;
    origin2[1] = ny-1;

    std::vector<T> origin3(2,T());
    origin3[0] = 0;
    origin3[1] = wallThickness;

    std::vector<T> origin4(2,T());
    origin4[0] = nx-wallThickness-1;
    origin4[1] = wallThickness;

    std::vector<T> extend1(2,T());
    extend1[0] = nx;
    extend1[1] = wallThickness;

    std::vector<T> extend2(2,T());
    extend2[0] = wallThickness;
    extend2[1] = ny-wallThickness;

    std::vector<T>extend3(2,T());
    extend3[0] = nx;
    extend3[1] = 1;

    Vector<T,2> centerDroplet(centerDropletX,centerDropletY);

    IndicatorCuboid2D<T> bottom(extend1, origin1);
    IndicatorCuboid2D<T> top(extend3, origin2);
    IndicatorCuboid2D<T> leftwall(extend2, origin3);
    IndicatorCuboid2D<T> rightwall(extend2, origin4);
    IndicatorCircle2D<T> droplet(centerDroplet, radiusDroplet);

    superGeometry.rename(1,2,droplet);
    superGeometry.rename(1,3,bottom);
    superGeometry.rename(1,3,leftwall);
    superGeometry.rename(1,3,rightwall);
    superGeometry.rename(1,4,top);

    /// Removes all not needed boundary voxels inside the surface
    superGeometry.innerClean();
    superGeometry.checkForErrors();

    superGeometry.print();

    clout << “Prepare Geometry … OK” << std::endl;
    }

    void prepareLattice(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    Dynamics<T, DESCRIPTOR>& bulkDynamics1,
    Dynamics<T, DESCRIPTOR>& bulkDynamics2,
    Dynamics<T, DESCRIPTOR>& bounceBackRho0,
    Dynamics<T, DESCRIPTOR>& bounceBackRho1,
    [b] sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeOne,
    sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeOne,
    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeTwo,
    sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeTwo,[/b]
    SuperGeometry2D<T>& superGeometry)
    {

    OstreamManager clout(std::cout,”prepareLattice”);
    clout << “Prepare Lattice …” << std::endl;

    [b] const T omega1 = 1.0;
    const T omega2 = 1.0;[/b]

    /// define lattice Dynamics
    sLatticeOne.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
    sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
    sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
    sLatticeOne.defineDynamics(superGeometry, 3, &bounceBackRho1);
    sLatticeOne.defineDynamics(superGeometry, 4, &bulkDynamics1);

    [b] sBoundaryConditionLatticeOne.addPressureBoundary(superGeometry, 4, omega1);[/b]

    sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
    sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
    sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
    sLatticeTwo.defineDynamics(superGeometry, 3, &bounceBackRho1);
    sLatticeTwo.defineDynamics(superGeometry, 4, &bulkDynamics2);
    [b]
    sBoundaryConditionLatticeTwo.addPressureBoundary(superGeometry, 4, omega2);[/b]

    clout << “Prepare Lattice … OK” << std::endl;
    }

    void setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    T force, int iT, SuperGeometry2D<T>& superGeometry) {

    if(iT==0) {

    AnalyticalConst2D<T,T> noise(4.e-2);
    std::vector<T> v(2,T());
    AnalyticalConst2D<T,T> zeroV(v);
    AnalyticalConst2D<T,T> zero(0.);
    AnalyticalLinear2D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.98+force*ny*DESCRIPTOR<T>::invCs2);
    AnalyticalConst2D<T,T> onePlus(0.98+force*ny/2.*DESCRIPTOR<T>::invCs2);
    AnalyticalRandom2D<T,T> random;
    AnalyticalIdentity2D<T,T> randomOne(random*noise+one);
    AnalyticalIdentity2D<T,T> randomPlus(random*noise+onePlus);
    std::vector<T> F(2,T());
    F[1] = -force;
    AnalyticalConst2D<T,T> f(F);

    /// for each material set the defineRhou and the Equilibrium

    sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 1,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 2,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeOne.defineRhoU(superGeometry, 4, zero, zeroV);
    sLatticeOne.iniEquilibrium(superGeometry, 4, zero, zeroV);
    sLatticeOne.defineExternalField(superGeometry, 4,
    DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
    DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

    sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);

    sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);

    sLatticeTwo.defineRhoU(superGeometry, 4, randomPlus, zeroV);
    sLatticeTwo.iniEquilibrium(superGeometry, 4, randomPlus, zeroV);

    /// Make the lattice ready for simulation
    sLatticeOne.initialize();
    sLatticeTwo.initialize();
    }
    }

    void getResults(SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
    SuperLattice2D<T, DESCRIPTOR>& sLatticeOne, int iT,
    SuperGeometry2D<T>& superGeometry, Timer<T>& timer) {

    OstreamManager clout(std::cout,”getResults”);
    SuperVTKwriter2D<T> vtkWriter(“rayleighTaylor2dsLatticeOne”);

    const int vtkIter = 100;
    const int statIter = 10;

    if (iT==0) {
    /// Writes the geometry, cuboid no. and rank no. as vti file for visualization
    SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLatticeOne, superGeometry);
    SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLatticeOne);
    SuperLatticeRank2D<T, DESCRIPTOR> rank(sLatticeOne);
    vtkWriter.write(geometry);
    vtkWriter.write(cuboid);
    vtkWriter.write(rank);
    vtkWriter.createMasterFile();
    }

    // Get statistics
    if (iT%statIter==0 && iT > 0) {
    // Timer console output
    timer.update( iT );
    timer.printStep();

    clout << “averageRhoFluidOne=” << sLatticeOne.getStatistics().getAverageRho();
    clout << “; averageRhoFluidTwo=” << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
    }

    /// Writes the VTK files
    if (iT%vtkIter==0) {
    clout << “Writing VTK …” << std::endl;
    SuperLatticeVelocity2D<T, DESCRIPTOR> velocity(sLatticeOne);
    SuperLatticeDensity2D<T, DESCRIPTOR> density(sLatticeOne);
    vtkWriter.addFunctor( velocity );
    vtkWriter.addFunctor( density );
    vtkWriter.write(iT);

    BlockLatticeReduction2D<T, DESCRIPTOR> planeReduction(density);
    BlockGifWriter<T> gifWriter;
    gifWriter.write(planeReduction, iT, “density”);

    clout << “Writing VTK … OK” << std::endl;
    }
    }

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

    /// === 1st Step: Initialization ===

    olbInit(&argc, &argv);
    singleton::directories().setOutputDir(“./tmp/”);
    OstreamManager clout(std::cout,”main”);

    const T omega1 = 1.0;
    const T omega2 = 1.0;
    const T G = 3.;
    T force = -30./(T)ny/(T)ny;

    /// === 2nd Step: Prepare Geometry ===
    /// Instantiation of a cuboidGeometry with weights

    #ifdef PARALLEL_MODE_MPI
    CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, singleton::mpi().getSize());
    #else
    CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, 1);
    #endif

    cGeometry.setPeriodicity(false, false);

    HeuristicLoadBalancer<T> loadBalancer(cGeometry);

    SuperGeometry2D<T> superGeometry(cGeometry, loadBalancer, 2);

    prepareGeometry(superGeometry);

    /// === 3rd Step: Prepare Lattice ===

    SuperLattice2D<T, DESCRIPTOR> sLatticeOne(superGeometry);
    SuperLattice2D<T, DESCRIPTOR> sLatticeTwo(superGeometry);

    ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
    omega1, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );
    ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics2 (
    omega2, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );

    // A bounce-back node with fictitious density 1,
    // which is experienced by the partner fluid
    BounceBack<T, DESCRIPTOR> bounceBackRho1( 1. );
    // A bounce-back node with fictitious density 0,
    // which is experienced by the partner fluid
    BounceBack<T, DESCRIPTOR> bounceBackRho0( 0. );

    std::vector<T> rho0;
    rho0.push_back(1);
    rho0.push_back(1);
    PsiEqualsRho<T,T> interactionPotential;
    ShanChenForcedGenerator2D<T,DESCRIPTOR> coupling(G,rho0,interactionPotential);

    [b] // Pressure boundary
    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeOne(sLatticeOne);
    createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeOne);
    sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeOne(sLatticeOne);

    sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeTwo(sLatticeTwo);
    createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeTwo);
    sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeTwo(sLatticeTwo);[/b]

    sLatticeOne.addLatticeCoupling(superGeometry, 1, coupling, sLatticeTwo);
    sLatticeOne.addLatticeCoupling(superGeometry, 2, coupling, sLatticeTwo);

    prepareLattice(sLatticeOne, sLatticeTwo, bulkDynamics1, bulkDynamics2,
    bounceBackRho0, bounceBackRho1,[b] sBoundaryConditionLatticeOne, sOffBoundaryConditionLatticeOne, sBoundaryConditionLatticeTwo, sOffBoundaryConditionLatticeTwo,[/b] superGeometry);

    /// === 4th Step: Main Loop with Timer ===
    int iT = 0;
    clout << “starting simulation…” << endl;
    Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
    timer.start();

    for (iT=0; iT<maxIter; ++iT) {

    /// === 5th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues(sLatticeOne, sLatticeTwo, force, iT, superGeometry);

    /// === 6th Step: Collide and Stream Execution ===
    sLatticeOne.collideAndStream();
    sLatticeTwo.collideAndStream();

    sLatticeOne.communicate();
    sLatticeTwo.communicate();

    sLatticeOne.executeCoupling();

    /// === 7th Step: Computation and Output of the Results ===
    getResults(sLatticeTwo, sLatticeOne, iT, superGeometry, timer);
    }

    timer.stop();
    timer.printSummary();
    }

    #2457
    ivan
    Member

    .

Viewing 3 posts - 1 through 3 (of 3 total)
  • You must be logged in to reply to this topic.