Skip to content

Particle periodicity running with MPI.

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Particle periodicity running with MPI.

Viewing 11 posts - 1 through 11 (of 11 total)
  • Author
    Posts
  • #8518
    Rookie
    Participant

    Dear OpenLB developers,

    /* This file is part of the OpenLB library
    *
    * Copyright (C) 2014 Thomas Henn, Mathias J. Krause
    * E-mail contact: info@openlb.net
    * The most recent release of OpenLB can be downloaded at
    * <http://www.openlb.net/&gt;
    *
    * This program is free software; you can redistribute it and/or
    * modify it under the terms of the GNU General Public License
    * as published by the Free Software Foundation; either version 2
    * of the License, or (at your option) any later version.
    *
    * This program is distributed in the hope that it will be useful,
    * but WITHOUT ANY WARRANTY; without even the implied warranty of
    * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    * GNU General Public License for more details.
    *
    * You should have received a copy of the GNU General Public
    * License along with this program; if not, write to the Free
    * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
    * Boston, MA 02110-1301, USA.
    */

    #ifndef PERIODICBOUNDARY3D_H_
    #define PERIODICBOUNDARY3D_H_

    #include <math.h>
    #include <vector>

    namespace olb {

    template<typename T, template<typename U> class PARTICLETYPE>
    class ParticleSystem3D;

    /*
    * Particle boundary based on a cube around the area with material number 1.
    * Only applicable to rectangles since if a particle leaves the area with
    * material number 1 it is moved to the opposing side of the area by
    * newPosition = oldPosition +/- extend(MaterialNumber=1).
    **/

    template<typename T, template<typename U> class PARTICLETYPE>
    class PeriodicBoundary3D : public Boundary3D<T, PARTICLETYPE> {
    public:
    PeriodicBoundary3D(SuperGeometry<T,3>& sg, bool x,
    bool y, bool z);
    PeriodicBoundary3D(PeriodicBoundary3D<T, PARTICLETYPE>& f);
    virtual ~PeriodicBoundary3D() { };
    virtual void applyBoundary(typename std::deque<PARTICLETYPE<T> >::iterator& p, ParticleSystem3D<T, PARTICLETYPE>& psSys);
    /// Returns number of particles that moved through the periodic boundary
    /// Order: x+, x-, y+, y-, z+, z-
    unsigned int* getJumper();
    private:
    //cube extents with origin (0,0,0)
    // std::vector<T> _minPhys, _maxPhys, _extend;
    olb::Vector<T, 3> _minPhys, _maxPhys, _extend;
    bool _x, _y, _z;
    unsigned int _jumper[6];
    CuboidGeometry3D<T>& _cuboidGeometry;
    T _overlap;
    };

    template<typename T, template<typename U> class PARTICLETYPE>
    PeriodicBoundary3D<T, PARTICLETYPE>::PeriodicBoundary3D(
    SuperGeometry<T,3>& sg, bool x, bool y, bool z) : Boundary3D<T, PARTICLETYPE>(),
    _minPhys(sg.getStatistics().getMinPhysR(1)),
    _maxPhys(sg.getStatistics().getMaxPhysR(1)),
    _extend(_maxPhys – _minPhys),
    _x(x),
    _y(y),
    _z(z), _cuboidGeometry(sg.getCuboidGeometry())
    {
    _minPhys = sg.getStatistics().getMinPhysR(1);
    _maxPhys = sg.getStatistics().getMaxPhysR(1);
    _extend[0] = _maxPhys[0] – _minPhys[0];
    _extend[1] = _maxPhys[1] – _minPhys[1];
    _extend[2] = _maxPhys[2] – _minPhys[2];
    for (int i=0; i<6; ++i) {
    _jumper[i] = 0;
    }
    _overlap = sg.getOverlap();
    }

    template<typename T, template<typename U> class PARTICLETYPE>
    void PeriodicBoundary3D<T, PARTICLETYPE>::applyBoundary(
    typename std::deque<PARTICLETYPE<T> >::iterator& p,
    ParticleSystem3D<T, PARTICLETYPE>& psSys)
    {
    if (_x) {
    if (p->getPos()[0] > _maxPhys[0]) {
    p->getPos()[0] -= _extend[0];
    ++_jumper[0];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    else if (p->getPos()[0] < _minPhys[0]) {
    p->getPos()[0] += _extend[0];
    ++_jumper[1];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    }
    if (_y) {
    if (p->getPos()[1] > _maxPhys[1]) {
    p->getPos()[1] -= _extend[1];
    ++_jumper[2];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    else if (p->getPos()[1] < _minPhys[1]) {
    p->getPos()[1] += _extend[1];
    ++_jumper[3];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    }
    if (_z) {
    if (p->getPos()[2] > _maxPhys[2]) {
    p->getPos()[2] -= _extend[2];
    ++_jumper[4];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    else if (p->getPos()[2] < _minPhys[2]) {
    p->getPos()[2] += _extend[2];
    ++_jumper[5];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);
    }
    }
    }

    template<typename T, template<typename U> class PARTICLETYPE>
    unsigned int* PeriodicBoundary3D<T, PARTICLETYPE>::getJumper()
    {
    return _jumper;
    }

    }

    #endif

    This is the code for particle periodic boundaries that I want to use. This code can only be used in serial mode, and when run in parallel, it encounters issues due to particles crossing block boundaries defined by MPI core count. To illustrate the problem, I have set up two particles, A and B, starting at positions A1 and B1 respectively. When they move to positions A2 and B2, the simulation terminates because particle A cannot enter the next block. I understand that this is because the code lacks relevant MPI settings. I have also found your documentation on the new settings for particle periodic boundaries (olb-1.6r0/src/particles/communication/relocation.h), which explains particle serialization and MPI non-blocking communication settings. The crux of the issue is that my code lacks MPI communication-related settings. I hope you can provide me with some guidance on how to modify it. I have tried adding MPI communication settings myself, but without success. I am very grateful for your assistance.

    https://postimg.cc/F7vmCHDM

    Best regards,
    Rookie

    #8570
    jan
    Participant

    Dear Rookie,

    are you using the subgrid3DLegacyFramework in your setup? Unfortunately, I cannot provide support in that context, so please correct me if I’m wrong.

    Here a few thoughts:

    I don’t see why the explained error should be caused by missing communications (I don’t think you have to add any communication, without the periodic boundary the particles move from one to the other block without problem, don’t they?). Because, if I understand it correctly, then the particles are duplicated when they pass through a block boundary in the middle of the domain. The periodic boundary shouldn’t be called there, should it?

    I could imagine that the hardcoded //cube extents with origin (0,0,0) doesn’t fit in your setup. Apparently, this part of the legacy code assumes that the simulation domain’s origin is at (0,0,0). Perhaps you should try to change your simulation setup to fit that assumption, if this doesn’t align with your setup.

    In general, consider to check these values (inside the applyBoundary method) for correctness, that is, are the values would they should be?

    _minPhys = sg.getStatistics().getMinPhysR(1);
    _maxPhys = sg.getStatistics().getMaxPhysR(1);
    _extend[0] = _maxPhys[0] – _minPhys[0];
    _extend[1] = _maxPhys[1] – _minPhys[1];
    _extend[2] = _maxPhys[2] – _minPhys[2];

    However, I’m aware that this doesn’t explain why it’s working in sequential mode, but not in parallel. Are you sure that it worked in sequential? Do the above values differ between sequential and parallel mode? I believe that here the positions are merely updated and that the communication should take place somewhere else.

    You could also try to change the places where the particle position is updated (in case some change to CuboidGeometry broke the logic, for example:

    p->getPos()[0] -= _extend[0];
    ++_jumper[0];
    int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    p->setCuboid(C);

    Consider following movePositionToEnd and movePositionToStart of the current code: https://www.openlb.net/DoxyGen/html/d9/d46/namespaceolb_1_1particles_1_1communication.html#a1647f5b5b28d226caa8ae64e54f0320b
    (min and max area already read from the SuperGeometry at the top and position refers to the current position)

    Best regards,
    Jan

    • This reply was modified 2 weeks, 5 days ago by jan.
    #8572
    jan
    Participant

    ~You could also try to change the places where the particle position is updated~
    You could also try to change the code where the particle position is updated

    #8576
    Rookie
    Participant

    Dear Jan,

    Thank you for answering my question. I am indeed using the subgrid3DLegacyFramework in my code. Although these codes are not easy to use, they do contain the functionality I need. Why aren’t you using this part of the code anymore? If this makes you uncomfortable, I won’t bother you with this issue again. I want to thank you for telling me not to add additional communication. From the results of inserting two particles, it seems that particles cannot move from one block to another. The problem may be that when particles pass through the boundary of the domain, they are not copied. Particle periodic boundaries should indeed not be called there. The question about whether the origin (0,0,0) for fluid and particle periodicity is set consistently. I set the origin of the fluid simulation domain to (0,0,0), but because of the fluid periodicity setting, an extra layer is added outside the geometric region, so it becomes (-Δx, -Δy, -Δz). I output the extend and minPhys of the particle periodicity, and these values ​​do not change in both serial and parallel modes. I only set periodicity in the x and y directions, so the extend values ​​for these two directions are correct. I want to emphasize that periodicity can run successfully in serial mode, so could this be the reason for termination in parallel mode?

    [SuperGeometryStatistics3D] materialNumber=1; count=39520; minPhysR=(0,0,0.00266667); maxPhysR=(0.250667,0.0826667,0.0346667)
    [SuperGeometryStatistics3D] materialNumber=2; count=6080; minPhysR=(0,0,0); maxPhysR=(0.250667,0.0826667,0.0373333)
    [prepareGeometry] Dimension of the channel in meters: x = 0.250667 ; y = 0.0826667 ; z = 0.0373333
    [prepareGeometry] Prepare Geometry … OK
    [main] noOfCuboid=8
    [prepareLattice] Prepare Lattice …
    [setInitialConditions] Set initial conditions …
    [setInitialConditions] Set initial conditions … OK
    [prepareLattice] Prepare Lattice … OK
    [main] Prepare Particles …
    _extend[0]=0.250667
    _extend[1]=0.0826667
    _extend[2]=0.032
    (_minPhys[0],_minPhys[1],_minPhys[2])=(0,0,0.00266667)

    Best regards,
    Rookie

    #8591
    mathias
    Keymaster

    You should change to the new code. It is much better. It should not be a big deal to transfer the functionallites you need from the old one to the new code. Once, you have done that, please consider contributing it to the OpenLB community by adding it at https://gitlab.com/openlb/release .

    #8633
    jan
    Participant

    Dear Rookie,

    unfortunately, I don’t have the capacities to check the legacy code in detail and it’s highly likely that some functionalities don’t work as expected anymore due to (recent) changes. However, just a few more thoughts on this matter:

    How do you incorporate the periodic boundary? Like this?

    auto periodicBoundary = make_shared < PeriodicBoundary3D<T, PARTICLE>
                              > (superGeometry, 1, 1, 1);
    spSys.addBoundary(periodicBoundary);
    

    Did you set an overlap such as spSys.setOverlap(pRadius * 2);?

    Also, I believe you didn’t try the last point of my previous post yet, right?

    Best regards,
    Jan

    • This reply was modified 1 week, 5 days ago by jan.
    #8639
    Rookie
    Participant

    Dear Jan,

    For your last point of suggestion, since I don’t know how to call the applyPeriodicityToPosition part of the code and whether to add it directly to the main function or another location. The new periodicity settings you provided should be applied to the new particle system (SuperParticleSystem3D<T, PARTICLE> spSys(cuboidGeometry, loadBalancer, superGeometry)). However, the functionality of this part of the code cannot implement the reverse action force of particles on the fluid, so legacy code is used. The particle system I am using is (SuperParticleSystem3D<T, PARTICLE> supParticleSystem(superGeometry)), and regarding the setting of overlap, it may be that the two lattices are too large? However, this can still run in serial mode.

      auto materialperiodicBoundary = std::make_shared
                                     < PeriodicBoundary3D<T, PARTICLE>
                                     > (superGeometry, true, true, false);
      supParticleSystem.addBoundary(materialperiodicBoundary); 
    
      supParticleSystem.setOverlap(2. * converter.getConversionFactorLength());
    
      auto dragModel = std::make_shared<SchillerNaumannDragModel<T, DESCRIPTOR, PARTICLE>>(converter);
      NaiveForwardCouplingModel<T, DESCRIPTOR, PARTICLE> forwardCoupling(converter, sLattice, superGeometry, dragModel);
      int overlap = 2;
      LocalBackCouplingModel<T, DESCRIPTOR, PARTICLE> backCoupling(converter, sLattice, superGeometry, overlap);
      int subSteps = 10;

    Best regards,
    Rookie

    #8688
    jan
    Participant

    Dear Rookie,

    how many blocks do you have in your system if you run the simulation sequential? Try using as many blocks as you do, when using the parallel system. Does the same problem occur in this case?

    It seems that you have a simplified case to work on this problem. Can you perhaps add an output of the particle position and cuboid at the end of applyBoundary? Does that align with your expectations? That could narrow down the possible causes.

    In general, I still don’t see how the periodic boundary could introduce the behavior shown in your image, because for me it seems that the periodic boundary shouldn’t be located at the intersection of these blocks. If it is, that might already be the problem that should be solved.
    Perhaps I misunderstand your setup. Could you provide some details or perhaps a sketch of the basic domain? Where are the periodic boundaries? How does the particle move in the images above, just from left to right?

    Is the fluid also periodic?

    Best regards,
    Jan

    • This reply was modified 5 days, 5 hours ago by jan.
    #8697
    Rookie
    Participant

    Dear Jan,

    Thank you for your patience with this issue. Below is my code. To speed up the simulation for you to test, I’ve shortened the simulation time. Please note that you must change the periodic boundaries to the code I posted earlier in order for it to compile successfully. If you enable parallel mode, you can join me in discovering issues together. Actually, regarding this code, I found that my “geometry2” should also contain particles, and particle collisions should occur in a layer of cells outside “geometry2”. However, it seems that the bidirectional coupling setting cannot be applied to two regions simultaneously.

    
    #include "olb3D.h"
    #include "olb3D.hh" // include full template code
    
    using namespace olb;
    using namespace olb::descriptors;
    using namespace olb::util;
    using namespace olb::graphics;
    
    typedef double T;
    typedef WallFunctionForcedD3Q19Descriptor DESCRIPTOR;
    #define PARTICLE Particle3D
    
    #ifndef M_PI
    #define M_PI 3.14159265358979323846
    #endif
    
    // Mathmatical constants
    const T pi = util::acos(-1);
    
    // Parameters for the simulation setup
    const int N = 59;
    const T physRefL = 0.02;       // half channel height in meters
    const T lx = 2. * pi * physRefL; // streamwise length in meters
    const T ly = 2. * physRefL;      // spanwise length in meters
    const T lz = 2. * physRefL;      // wall-normal length in meters
    const T radius = 2.5e-5;        // particles radius
    const T partRho = 2500.;         // particles density
    
    // Choose friction reynolds number ReTau
    // #define Case_ReTau_1000
    // #define Case_ReTau_2000
    #define Case_ReTau_180
    
    // Wallfunction parameters
    const T latticeWallDistance = 0.5; // lattice distance to boundary
    const int rhoMethod = 2;           // method for density reconstruction
    // 0: Zou-He
    // 1: extrapolation
    // 2: constant
    const int fneqMethod = 0; // method for fneq reconstruction
    // 0: regularized NEBB (Latt)
    // 1: extrapolation NEQ (Guo Zhaoli)
    // 2: regularized second order finite Differnce
    // 3: equilibrium scheme
    const int wallProfile = 0; // wallfunction profile
    // 0: Musker profile
    // 1: power law profile
    
    // Reynolds number based on the friction velocity
    #if defined(Case_ReTau_1000)
    T ReTau = 1000.512;
    #elif defined(Case_ReTau_2000)
    T ReTau = 1999.756;
    #elif defined(Case_ReTau_180)
    T ReTau = 180;
    #endif
    // Characteristic physical kinematic viscosity from DNS Data
    // http://turbulence.ices.utexas.edu/channel2015/data/LM_Channel_1000_mean_prof.dat
    #if defined(Case_ReTau_1000)
    T charPhysNu = 5. / 100000.;
    #elif defined(Case_ReTau_2000)
    T charPhysNu = 2.3 / 100000.;
    #elif defined(Case_ReTau_180)
    T charPhysNu = 1.5 / 100000.;
    #endif
    
    // number of forcing updates over simulation time
    #if defined(Case_ReTau_1000)
    T fluxUpdates = 2000;
    #elif defined(Case_ReTau_2000)
    T fluxUpdates = 4000;
    #elif defined(Case_ReTau_180)
    T fluxUpdates = 2000;
    #endif
    
    // physical simulated length adapted for lattice distance to boundary in meters
    const T adaptedPhysSimulatedLength = 2 * physRefL - 2 * ((0.04 / T(N + 2 * latticeWallDistance)) * latticeWallDistance);
    const T Re = 6594;
    // Characteristic physical mean bulk velocity from Dean correlations in meters - Malaspinas and Sagaut (2014)
    const T charPhysU = Re * charPhysNu / (2. * physRefL);
    // Time of the simulation in seconds
    const T charPhysT = physRefL / (ReTau * charPhysNu / physRefL);
    
    const T physConvergeTime = 12. * charPhysT;                   // time until until statistics sampling in seconds
    const T physStatisticsTime = 4. * charPhysT;                  // statistics sampling time in seconds
    const T particleMaxPhysT = 1. * charPhysT;
    const T singleMaxPhysT = physConvergeTime + physStatisticsTime;
    const T fluidMaxPhysT = physConvergeTime + physStatisticsTime + particleMaxPhysT; // max. simulation time in seconds
    const T statisticsSave = 1. / 25.; // time between statistics samples in seconds
    const int noOfParticles = 30000;
    const T checkstatistics = (T)fluidMaxPhysT / 200.;
    
    // seed the rng with time if SEED_WITH_TIME is set, otherwise just use a fixed seed.
    #if defined SEED_WITH_TIME
    #include <chrono>
    auto seed = std::chrono::system_clock().now().time_since_epoch().count();
    std::default_random_engine generator(seed);
    #else
    std::default_random_engine generator(0x1337533DAAAAAAAA);
    #endif
    
    // Compute mean lattice velocity from musker wallfunction
    T computeLatticeVelocity()
    {
      T Ma_max = 0.1;
      T c_s = 1 / util::sqrt(3.0);
      T latticeUMax = Ma_max * c_s;
      Musker<T, T> musker_tmp(charPhysNu, physRefL, 1.205);
      T charPhysU_tau = ReTau * charPhysNu / physRefL;
      T tau_w[1];
      tau_w[0] = util::pow(charPhysU_tau, 2.);
      T charPhysUMax[1];
      musker_tmp(charPhysUMax, tau_w);
      T latticeU = charPhysU * latticeUMax / charPhysUMax[0];
      return latticeU;
    }
    
    template <typename T, typename S>
    class Channel3D : public AnalyticalF3D<T, S>
    {
    
    protected:
      T turbulenceIntensity;
      T maxVelocity;
      T distanceToWall;
      T obst_z;
      T obst_r;
      T a;
      T b;
    
    public:
      Channel3D(UnitConverter<T, DESCRIPTOR> const &converter, T frac) : AnalyticalF3D<T, S>(3)
      {
        turbulenceIntensity = 0.05;
        distanceToWall = -converter.getPhysDeltaX() / 2.;
        maxVelocity = converter.getLatticeVelocity(converter.getCharPhysVelocity() * (8. / 7.)); // Centerline Velocity
        obst_z = physRefL + distanceToWall;
        obst_r = physRefL;
        a = -1.;
        b = 1.;
      };
    
      bool operator()(T output[], const S input[])
      {
        std::uniform_real_distribution<T> distribution(a, b);
        T nRandom1 = distribution(generator);
        T nRandom2 = distribution(generator);
        T nRandom3 = distribution(generator);
    
        T u_calc = maxVelocity * util::pow(((obst_r - util::abs(input[2] - obst_z)) / obst_r), 1. / 7.);
    
        output[0] = turbulenceIntensity * nRandom1 * maxVelocity + u_calc;
        output[1] = turbulenceIntensity * nRandom2 * maxVelocity;
        output[2] = turbulenceIntensity * nRandom3 * maxVelocity;
    
        return true;
      };
    };
    
    template <typename T, typename S>
    class TrackedForcing3D : public AnalyticalF3D<T, S>
    {
    
    protected:
      T um;
      T utau;
      T h2;
      T aveVelocity;
    
    public:
      TrackedForcing3D(UnitConverter<T, DESCRIPTOR> const &converter, int ReTau) : AnalyticalF3D<T, S>(3)
      {
        um = converter.getCharPhysVelocity();
        utau = ReTau * converter.getPhysViscosity() / (converter.getCharPhysLength() / 2.);
        h2 = converter.getCharPhysLength() / 2.;
        aveVelocity = um;
      };
    
      void updateAveVelocity(T newVel)
      {
        aveVelocity = newVel;
      }
    
      bool operator()(T output[], const S input[])
      {
        output[0] = util::pow(utau, 2) / h2 + (um - aveVelocity) * um / h2;
        output[1] = 0;
        output[2] = 0;
    
        return true;
      };
    };
    
    void prepareGeometry(SuperGeometry<T, 3> &superGeometry, IndicatorF3D<T> &indicator,
                         UnitConverter<T, DESCRIPTOR> const &converter)
    {
      OstreamManager clout(std::cout, "prepareGeometry");
      clout << "Prepare Geometry ..." << std::endl;
    
      superGeometry.rename(0, 2, indicator);
      superGeometry.rename(2, 1, {0, 0, 1});
      superGeometry.clean();
      superGeometry.innerClean();
      superGeometry.checkForErrors();
    
      superGeometry.print();
    
      olb::Vector<T, 3> PhyMax = superGeometry.getStatistics().getMaxPhysR(2);
      olb::Vector<T, 3> PhyMin = superGeometry.getStatistics().getMinPhysR(2);
      clout << "Dimension of the channel in meters: x = " << PhyMax[0] - PhyMin[0];
      clout << " ; y = " << PhyMax[1] - PhyMin[1];
      clout << " ; z = " << PhyMax[2] - PhyMin[2] << std::endl;
    
      clout << "Prepare Geometry ... OK" << std::endl;
    }
    
    // set up initial conditions
    void setInitialConditions(SuperLattice<T, DESCRIPTOR> &sLattice,
                              UnitConverter<T, DESCRIPTOR> const &converter,
                              SuperGeometry<T, 3> &superGeometry,
                              AnalyticalScaled3D<T, T> &forceSolScaled,
                              TrackedForcing3D<T, T> &forceSol)
    {
    
      OstreamManager clout(std::cout, "setInitialConditions");
      clout << "Set initial conditions ..." << std::endl;
    
      AnalyticalConst3D<T, T> rho(1.205);
      Channel3D<T, T> uSol(converter, 1.);
    
      sLattice.defineRhoU(superGeometry, 1, rho, uSol);
      sLattice.iniEquilibrium(superGeometry, 1, rho, uSol);
    
      sLattice.defineRhoU(superGeometry, 2, rho, uSol);
      sLattice.iniEquilibrium(superGeometry, 2, rho, uSol);
    
      AnalyticalConst3D<T, T> TauEff(1. / converter.getLatticeRelaxationFrequency());
    
      sLattice.defineField<TAU_EFF>(superGeometry, 1, TauEff);
      sLattice.defineField<TAU_EFF>(superGeometry, 2, TauEff);
    
      // Force Initialization
      forceSol.updateAveVelocity(converter.getCharPhysVelocity()); // New average velocity
    
      // Initialize force
      sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
      sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
      // Tau_w Initialization
      T tau_w_guess = 0.0; // Wall shear stress in phys units
      AnalyticalConst3D<T, T> tau_w_ini(tau_w_guess);
      AnalyticalScaled3D<T, T> tau_w_ini_scaled(tau_w_ini, 1. / (converter.getConversionFactorForce() * util::pow(converter.getConversionFactorLength(), 2.)));
    
      sLattice.defineField<TAU_W>(superGeometry, 1, tau_w_ini_scaled);
      sLattice.defineField<TAU_W>(superGeometry, 2, tau_w_ini_scaled);
    
      clout << "Set initial conditions ... OK" << std::endl;
    }
    
    // Set up the geometry of the simulation
    void prepareLattice(SuperLattice<T, DESCRIPTOR> &sLattice,
                        UnitConverter<T, DESCRIPTOR> const &converter,
                        SuperGeometry<T, 3> &superGeometry,
                        AnalyticalScaled3D<T, T> &forceSolScaled,
                        TrackedForcing3D<T, T> &forceSol,
                        wallFunctionParam<T> const &wallFunctionParam)
    {
    
      OstreamManager clout(std::cout, "prepareLattice");
      clout << "Prepare Lattice ..." << std::endl;
    
      /// Material=1 -->bulk dynamics
      sLattice.defineDynamics<SmagorinskyForcedBGKdynamics>(superGeometry, 1);
      /// Material = 2 --> boundary node + wallfunction
      sLattice.defineDynamics<ExternalTauEffLESForcedBGKdynamics>(superGeometry, 2);
      setWallFunctionBoundary<T, DESCRIPTOR>(sLattice, superGeometry, 2, converter, wallFunctionParam);
    
      /// === Set Initial Conditions == ///
      setInitialConditions(sLattice, converter, superGeometry, forceSolScaled, forceSol);
    
      sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());
      sLattice.setParameter<collision::LES::Smagorinsky>(0.1);
    
      // Make the lattice ready for simulation
      sLattice.initialize();
    
      clout << "Prepare Lattice ... OK" << std::endl;
    }
    
    /// Computes the pressure drop between the voxels before and after the cylinder
    bool getResults(SuperLattice<T, DESCRIPTOR> &sLattice,
                    UnitConverter<T, DESCRIPTOR> const &converter, size_t iT, int iTperiod,
                    SuperGeometry<T, 3> &superGeometry, Timer<double> &fluidTimer,
                    SuperParticleSystem3D<T, PARTICLE> &supParticleSystem,
                    T radii, T partRho, Timer<double> &particleTimer,
                    SuperParticleSysVtuWriter<T, PARTICLE> &supParticleWriter,
                    bool fluidExists, SuperLatticeTimeAveragedF3D<T> &sAveragedVel)
    {
      OstreamManager clout(std::cout, "getResults");
    
      std::list<int> materialslist;
      materialslist.push_back(1);
      materialslist.push_back(2);
      using BulkDynamics = ShearSmagorinskyBGKdynamics<T,DESCRIPTOR>;
      ParametersOfOperatorD<T,DESCRIPTOR,BulkDynamics> bulkDynamicsParams{};
    
      SuperVTMwriter3D<T> vtmWriter("channel3d");
      SuperVTMwriter3D<T> vtmWriterStartTime("startingTimechannel3d");
      SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
      SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
      SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
      vtmWriter.addFunctor(geometry);
      vtmWriter.addFunctor(velocity);
      vtmWriter.addFunctor(pressure);
    
      std::size_t singleMaxT = converter.getLatticeTime(singleMaxPhysT);
    
      if (iT == 0)
      {
        // Writes the geometry, cuboid no. and rank no. as vti file for visualization
        SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
        SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
        SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
    
        vtmWriter.write(geometry);
        vtmWriter.write(cuboid);
        vtmWriter.write(rank);
    
        vtmWriter.createMasterFile();
        vtmWriterStartTime.createMasterFile();
    
        // Print some output of the chosen simulation setup
        clout << "N=" << N << "; maxTimeSteps(fluid)="
              << converter.getLatticeTime(fluidMaxPhysT) << "; noOfCuboid="
              << superGeometry.getCuboidGeometry().getNc() << "; Re=" << Re
              << "; noOfParticles=" << noOfParticles << "; maxTimeSteps(particle)="
              << converter.getLatticeTime(particleMaxPhysT)
              << "; St=" << (2. * partRho * radius * radius * converter.getCharPhysVelocity()) / (9. * converter.getPhysViscosity() * converter.getPhysDensity() * converter.getCharPhysLength()) << std::endl;
      }
    
      // Writes output on the console for the fluid phase
      if (iT < converter.getLatticeTime(fluidMaxPhysT) && iT % iTperiod == 0)
      {
        // Timer console output
        fluidTimer.update(iT);
        fluidTimer.printStep(2);
        // Lattice statistics console output
        sLattice.getStatistics().print(iT, iT * converter.getPhysDeltaT());
        clout << "Max. physical velocity(m/s): " << converter.getPhysVelocity(sLattice.getStatistics().getMaxU()) << std::endl;
        clout << "Max u+:" << converter.getPhysVelocity(sLattice.getStatistics().getMaxU()) / (ReTau * charPhysNu / (converter.getCharPhysLength() / 2.)) << std::endl;
      }
    
      if (iT < converter.getLatticeTime(fluidMaxPhysT) && iT % converter.getLatticeTime(statisticsSave) == 0 && iT > converter.getLatticeTime(physConvergeTime))
      {
        // Add ensemble to temporal averaged velocity
        sLattice.communicate();
        sAveragedVel.addEnsemble();
      }
    
      if (iT < converter.getLatticeTime(singleMaxPhysT) && (iT % converter.getLatticeTime(singleMaxPhysT / 16) == 0 || iT == converter.getLatticeTime(fluidMaxPhysT) - 1))
      {
        // Writes the vtk files
        vtmWriter.write(iT);
      }
    
      // Writes output on the console for the fluid phase
      if (iT >= converter.getLatticeTime(singleMaxPhysT) &&
          (iT % (iTperiod) == 0  || iT == converter.getLatticeTime(singleMaxPhysT) || iT == converter.getLatticeTime(fluidMaxPhysT) - 1))
      {
        vtmWriter.write(iT);
        
        particleTimer.print(iT - singleMaxT);
    
        // console output number of particles at different material numbers mat
        supParticleSystem.print({1, 2});
    
        // only write .vtk-files after the fluid calculation is finished
        supParticleWriter.write(iT - singleMaxT);
    
        // true as long as certain amount of active particles
        if (supParticleSystem.globalNumOfActiveParticles() < 0.0001 * noOfParticles && iT > 0.9 * converter.getLatticeTime(fluidMaxPhysT))
        {
          return false;
        }
      }
      return true;
    }
    
    int main(int argc, char *argv[])
    {
    
      /// === 1st Step: Initialization ===
      olbInit(&argc, &argv);
      singleton::directories().setOutputDir("./tmp/");
      OstreamManager clout(std::cout, "main");
      // display messages from every single mpi process
      // clout.setMultiOutput(true);
    
      UnitConverterFromResolutionAndLatticeVelocity<T, DESCRIPTOR>  converter(
          int{N},                        // resolution: number of voxels per charPhysL
          (T)computeLatticeVelocity(),   // latticeU : mean lattice velocity
          (T)adaptedPhysSimulatedLength, // charPhysLength: reference length of simulation geometry
          (T)2.19,                        // charPhysVelocity: mean bulk velocity in __m / s__
          (T)charPhysNu,                 // physViscosity: physical kinematic viscosity in __m^2 / s__
          (T)1.205                         // physDensity: physical density in __kg / m^3__
      );
      converter.print();
      converter.write("channelpraticle");
    
      clout << "----------------------------------------------------------------------" << std::endl;
      clout << "Converge time(s): " << physConvergeTime << std::endl;
      clout << "Lattice converge time: " << converter.getLatticeTime(physConvergeTime) << std::endl;
      clout << "Max. Phys. simulation time(s): " << fluidMaxPhysT << std::endl;
      clout << "Max. Lattice simulation time: " << converter.getLatticeTime(fluidMaxPhysT) << std::endl;
      clout << "Frequency Statistics Save(Hz): " << 1. / statisticsSave << std::endl;
      clout << "Statistics save period(s): " << statisticsSave << std::endl;
      clout << "Lattice statistics save period: " << converter.getLatticeTime(statisticsSave) << std::endl;
      clout << "----------------------------------------------------------------------" << std::endl;
      clout << "Channel height(m): " << adaptedPhysSimulatedLength << std::endl;
      clout << "y+ value: " << (ReTau * converter.getPhysViscosity() / (physRefL)) * ((0.04 / T(N + 2 * latticeWallDistance)) * latticeWallDistance) / converter.getPhysViscosity() << std::endl;
      clout << "y+ value spacing: " << (ReTau * converter.getPhysViscosity() / (physRefL)) * (converter.getPhysDeltaX()) / converter.getPhysViscosity() << std::endl;
      clout << "----------------------------------------------------------------------" << std::endl;
      clout << "N=" << N << "; maxTimeSteps(fluid)="
            << converter.getLatticeTime(fluidMaxPhysT) << "; Re=" << Re
            << "; noOfParticles=" << noOfParticles << "; maxTimeSteps(particle)="
            << converter.getLatticeTime(particleMaxPhysT) << "; singleMaxPhysT="
            << converter.getLatticeTime(singleMaxPhysT)
            << "; St=" << (2. * partRho * radius * radius * converter.getCharPhysVelocity()) / (9. * converter.getPhysViscosity() * converter.getPhysDensity() * converter.getCharPhysLength()) << std::endl;
      clout << "----------------------------------------------------------------------" << std::endl;
      clout << "iTperiod=" << converter.getLatticeTime(checkstatistics) << std::endl;
      clout << "utau=" << ReTau * converter.getPhysViscosity() / (converter.getCharPhysLength() / 2.) << std::endl;
      clout << "----------------------------------------------------------------------" << std::endl;
    
    #ifdef PARALLEL_MODE_MPI
      const int noOfCuboids = singleton::mpi().getSize();
    #else
      const int noOfCuboids = 1;
    #endif
    
      Vector<T, 3> extend(lx, ly, adaptedPhysSimulatedLength);
      extend[2] += (1. / 8.) * converter.getPhysDeltaX();
    
      Vector<T, 3> origin(0., 0., 0.);
      IndicatorCuboid3D<T> cuboid(extend, origin);
    
      CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids);
    
      cuboidGeometry.setPeriodicity(true, true, false);
    
      HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
    
      SuperGeometry<T, 3> superGeometry(cuboidGeometry, loadBalancer, 3);
    
      prepareGeometry(superGeometry, cuboid, converter);
      clout << "noOfCuboid=" << superGeometry.getCuboidGeometry().getNc() << std::endl;
    
      /// === 3rd Step: Prepare Lattice ===
      SuperLattice<T, DESCRIPTOR> sLattice(superGeometry);
    
      // forcing of the channel
      TrackedForcing3D<T, T> forceSol(converter, ReTau);
      AnalyticalScaled3D<T, T> forceSolScaled(forceSol, 1. / (converter.getConversionFactorForce() / converter.getConversionFactorMass()));
    
      int input[3];
      T output[5];
    
      Vector<T, 3> normal(1, 0, 0);
      std::vector<int> normalvec{1, 0, 0};
      Vector<T, 3> center;
      for (int i = 0; i < 3; ++i)
      {
        center[i] = origin[i] + extend[i] / 2;
      }
      SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
      std::vector<int> materials;
      materials.push_back(1);
      materials.push_back(2);
    
      std::list<int> materialslist;
      materialslist.push_back(1);
      materialslist.push_back(2);
    
      wallFunctionParam<T> wallFunctionParam;
      wallFunctionParam.bodyForce = true;
      wallFunctionParam.wallProfile = wallProfile;
      wallFunctionParam.rhoMethod = rhoMethod;
      wallFunctionParam.fneqMethod = fneqMethod;
      wallFunctionParam.latticeWalldistance = latticeWallDistance;
      wallFunctionParam.vonKarman = 0.4;
      wallFunctionParam.curved = false;
    
      prepareLattice(sLattice, converter, superGeometry,
                     forceSolScaled, forceSol, wallFunctionParam);
    
      SuperPlaneIntegralFluxVelocity3D<T> velFlux(sLattice, converter, superGeometry, center, normal, materials,
                                                  BlockDataReductionMode::Discrete);
    
      // === 3.1 Step: Particles ===
      clout << "Prepare Particles ..." << std::endl;
    
      // SuperParticleSystems3D
      SuperParticleSystem3D<T, PARTICLE> supParticleSystem(superGeometry);
      // define which properties are to be written in output data
      SuperParticleSysVtuWriter<T, PARTICLE> supParticleWriter(supParticleSystem,
                                                               "particles", SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::velocity | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::mass | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::radius 
                                                               | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::active | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::force);
    
      SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR> getVel(sLattice, converter);
      T dynVisc = converter.getPhysViscosity() * converter.getPhysDensity();
      T physDensity = converter.getPhysDensity();
    
      auto schillerNaumannDragForce = std::make_shared<SchillerNaumannDragForce3D<T, PARTICLE, DESCRIPTOR>>(getVel, dynVisc, physDensity);
      supParticleSystem.addForce(schillerNaumannDragForce);
    
      std::vector<T> direction{1, 0, 0};
      const T g = 9.81;
      auto weightForce = std::make_shared<WeightForce3D<T, PARTICLE>>(direction, g);
      supParticleSystem.addForce(weightForce);
    
      T dT = converter.getConversionFactorTime();
      std::set<int> reflBMat = {2};
      auto materialreflectBoundary = std::make_shared<SimpleReflectBoundary3D<T, PARTICLE>>(dT, superGeometry, reflBMat);
      supParticleSystem.addBoundary(materialreflectBoundary);
    
      auto materialperiodicBoundary = std::make_shared
                                     < PeriodicBoundary3D<T, PARTICLE>
                                     > (superGeometry, true, true, false);
      supParticleSystem.addBoundary(materialperiodicBoundary); 
    
      supParticleSystem.setOverlap(2. * converter.getConversionFactorLength());
     
      auto dragModel = std::make_shared<SchillerNaumannDragModel<T, DESCRIPTOR, PARTICLE>>(converter);
      NaiveForwardCouplingModel<T, DESCRIPTOR, PARTICLE> forwardCoupling(converter, sLattice, superGeometry, dragModel);
      int overlap = 2;
      LocalBackCouplingModel<T, DESCRIPTOR, PARTICLE> backCoupling(converter, sLattice, superGeometry, overlap);
      int subSteps = 10;
    
      // particles generation at inlet3
      Vector<T, 3> extendP = {lx, ly, adaptedPhysSimulatedLength};
      Vector<T, 3> originP = {0, 0, 0};
      IndicatorCuboid3D<T> inletCuboid(extendP, originP);
      supParticleSystem.addParticle(inletCuboid, 4. / 3. * M_PI * util::pow(radius, 3) * partRho, radius, noOfParticles);
    
      clout << "Prepare Particles ... OK" << std::endl;
    
      /// === 5th Step: Definition of turbulent Statistics Objects ===
    
      SuperLatticePhysVelocity3D<T, DESCRIPTOR> sVel(sLattice, converter);
      SuperLatticeTimeAveragedF3D<T> sAveragedVel(sVel);
      SuperLatticePhysPressure3D<T, DESCRIPTOR> sPre(sLattice, converter);
      
      /// === 4th Step: Main Loop with Timer ===
      Timer<double> fluidTimer(converter.getLatticeTime(fluidMaxPhysT), superGeometry.getStatistics().getNvoxel());
    
      Timer<double> particleTimer(converter.getLatticeTime(particleMaxPhysT),
                                  noOfParticles);
    
      fluidTimer.start();
    
      std::size_t iT = 0;
    
      int iTperiod = converter.getLatticeTime(checkstatistics);
    
      bool fluidExists = true;
    
      // checks whether there is already data of the fluid from an earlier calculation
      if (!(sLattice.load("fluidSolution")))
      {
    
        fluidExists = false;
    
        for (; iT <= converter.getLatticeTime(singleMaxPhysT); ++iT)
        {
          if (iT % converter.getLatticeTime(fluidMaxPhysT / fluxUpdates) == 0 || iT == 0)
          {
    
            velFlux(output, input);
            T flux = output[0];
            T area = output[1];
            forceSol.updateAveVelocity(flux / area);
            sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
            sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
          }
          
          /// === 6th Step: Computation and Output of the Results ===
          getResults(sLattice, converter, iT, iTperiod, superGeometry, fluidTimer,
                     supParticleSystem, radius, partRho, particleTimer,
                     supParticleWriter, fluidExists, sAveragedVel);
    
          
    
          /// === 7th Step: Collide and Stream Execution ===
          sLattice.collideAndStream();
        }
          sLattice.save("fluidSolution");
      }
      else
      {
          iT = converter.getLatticeTime(singleMaxPhysT);
          getResults(sLattice, converter, iT,
                     iTperiod, superGeometry, fluidTimer,
                     supParticleSystem, radius, partRho, particleTimer,
                     supParticleWriter, fluidExists, sAveragedVel);
      }   
          // when the fluid simulation time reaches singleMaxPhysT, the particle simulation begins
          supParticleSystem.setVelToFluidVel( getVel );
          particleTimer.start();
    
          for ( ; iT <= converter.getLatticeTime( fluidMaxPhysT ); ++iT ) {
            
            if (iT % converter.getLatticeTime(fluidMaxPhysT / fluxUpdates) == 0)
          {
    
            velFlux(output, input);
            T flux = output[0];
            T area = output[1];
            forceSol.updateAveVelocity(flux / area);
            sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
            sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
          }
            // particles simulation starts after run up time is over
            // supParticleSystem.simulate( converter.getConversionFactorTime());
            supParticleSystem.simulateWithTwoWayCoupling_Mathias ( dT, forwardCoupling, backCoupling, 1, subSteps, true);
    
            if ( !getResults( sLattice, converter, iT,
                              iTperiod, superGeometry, fluidTimer,
                              supParticleSystem, radius, partRho, particleTimer,
                              supParticleWriter, fluidExists, sAveragedVel) ) 
              {
                break;
              }
             
          /// === 7th Step: Collide and Stream Execution ===
          sLattice.collideAndStream();
    
      }
    
       fluidTimer.stop();
       fluidTimer.printSummary();
       particleTimer.stop();
       particleTimer.printSummary();
    }
    
    

    Best regards,
    Rookie

    • This reply was modified 4 days, 9 hours ago by Rookie.
    #8705
    jan
    Participant

    Dear Rookie,

    Unfortunately, I cannot run the provided code because it seems to be missing some parts. However, I suggest that you use a minimal working example where you experience the problem anyway, so that it’s easier to identify the actual problem. If it’s in the periodic boundaries, then a very simple setup with a particle moving at constant speed through the periodic boundary should suffice, right?

    Best regards,
    Jan

    #8715
    Rookie
    Participant

    Dear Jan,

    First, you need to run it under version 1.6. Then, can you tell me the errors you encountered during compilation or any issues that occurred during runtime? I’ll see if I can solve them. If not, let’s just give up on this problem.

    Best regards,
    Rookie

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