[Rookie]

Dear Christoph,

Thanks for your suggestion, I have solved this problem, because this function（VerletParticleDynamicsVelocityWallReflection） only accepts one parameter, so defining two walls is not successful, so I changed the size of the cuboid.

Best regards,
Rookie

[Rookie]

Dear Jan,

I want the resolved particles to collide with the walls in the channel, so I need solid boundaries at both the top and bottom. My fluid and particles are set to be periodic in the x and y directions. However, if I use the code mentioned earlier, it indeed causes particle rebound in the x and y directions due to walls surrounding all sides. The problem now is that I don’t know how to properly set up the walls.

I only found “SuperIndicatorMaterial,” but it cannot be directly used to define the walls due to material mismatches. As you mentioned, “IndicatorCuboid” seems to define only one side of the wall, but I need walls on both the top and bottom. Regarding your suggestion to provide sufficient thickness, is it reflected in “5 * converter.getPhysDeltaX()”? However, based on my tests, particles rebound off the surrounding walls, and changing this parameter does not affect the simulation results.

const unsigned latticeMaterial = 2; //Material number of wall
const unsigned contactMaterial = 0; //Material identifier (only relevant for contact model)
SolidBoundary<T,3> wall( std::make_unique<IndicInverse<T, DESCRIPTOR::d>>(
cuboid, cuboid.getMin() – 5 * converter.getPhysDeltaX(),
cuboid.getMax() + 5 * converter.getPhysDeltaX()),
latticeMaterial, contactMaterial );

particleSystem.defineDynamics<
VerletParticleDynamicsVelocityWallReflection<T,PARTICLETYPE>>(wall);

https://postimg.cc/G9qCX7Bb

Best regards,
Rookie

[Rookie]

Dear Jan,

I’ve noticed an interesting phenomenon: the reduction in the number of particles only occurs on servers with AMD CPUs. However, my personal computer and another server both have Intel CPUs, and the number of particles doesn’t decrease. I’ve checked the code, and it’s the same across all systems. Could this be the cause of the issue?

Best regards,
Rookie

[Rookie]

The following picture sets the simulation results of two particles in 8 threads, and it is obvious that the periodic boundary of the particles has been achieved.
https://postimg.cc/PpdV8VX5

[Rookie]

Dear Jan,

Yes, deleting this sentence will allow the simulation to continue. If not, an error will be reported, and the particle will move around the periodic boundary after modifying the particle’s position coordinates. However, I found a new problem. If I increase the number of particles and the number of CPU cores, it will cause new problems in MPI communication, such as the number of particles will decrease a little, or the following error [ubuntu:14564] [[15031,0],0] ORTE_ERROR_LOG: Data unpack would read past end of buffer in file util/show_help.c at line 507,is this related to the fact that I no longer transfer these particles to the neighboring block? At present, I am still learning MPI related knowledge, if it is solved, I will continue to share it.

Best regards,
Rookie

[Rookie]

Dear Jan,

I am very glad that I have fixed the code for the periodic boundary. However, I am not sure if I fixed it correctly. Initially, I thought it was an issue with MPI non-blocking communication, so I searched for information and tried using MPI_Barrier(MPI_COMM_WORLD). I found that the program could run correctly when the number of particles did not exceed the simulation domain, but it would deadlock when particles crossed the boundary. Therefore, I tried commenting out p->setCuboid(newCuboid) because particle allocation was no longer necessary. This confirmed that I correctly implemented the periodic boundary for particles; it prints that particles move from the current rank to the opposite side of the new cuboid. Could you help me confirm if I made the correct modifications? You can quickly complete the simulation by slightly reducing the simulation time and the number of particles.

#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)
  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)
{
  bool crossed_boundary = false;
  
  if (_x) {
    if (p->getPos()[0] > _maxPhys[0]) {
      p->getPos()[0] -= _extend[0];
      ++_jumper[0];
      crossed_boundary = true;
      std::cout << "Particle crossed x+ boundary." << std::endl;
    }
    else if (p->getPos()[0] < _minPhys[0]) {
      p->getPos()[0] += _extend[0];
      ++_jumper[1];
      crossed_boundary = true;
      std::cout << "Particle crossed x- boundary." << std::endl;
    }
  }
  if (_y) {
    if (p->getPos()[1] > _maxPhys[1]) {
      p->getPos()[1] -= _extend[1];
      ++_jumper[2];
      crossed_boundary = true;
      std::cout << "Particle crossed y+ boundary." << std::endl;
    }
    else if (p->getPos()[1] < _minPhys[1]) {
      p->getPos()[1] += _extend[1];
      ++_jumper[3];
      crossed_boundary = true;
      std::cout << "Particle crossed y- boundary." << std::endl;
    }
  }
  if (_z) {
    if (p->getPos()[2] > _maxPhys[2]) {
      p->getPos()[2] -= _extend[2];
      ++_jumper[4];
      crossed_boundary = true;
      std::cout << "Particle crossed z+ boundary." << std::endl;
    }
    else if (p->getPos()[2] < _minPhys[2]) {
      p->getPos()[2] += _extend[2];
      ++_jumper[5];
      crossed_boundary = true;
      std::cout << "Particle crossed z- boundary." << std::endl;
    }
  }

  if (crossed_boundary) {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    int newCuboid = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    // p->setCuboid(newCuboid);

    std::cout << "Rank " << rank << ": Particle crossed boundary to new cuboid " << newCuboid << std::endl;

    // Ensure synchronization among all processes.
    // MPI_Barrier(MPI_COMM_WORLD);
  }
}

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

} // namespace olb

#endif

const T physConvergeTime = 0.1 * charPhysT;                   // time until until statistics sampling in seconds
const T physStatisticsTime = 0.1 * charPhysT;                  // statistics sampling time in seconds
const T particleMaxPhysT = 0.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 = 100;
const T checkstatistics = (T)fluidMaxPhysT / 200.;

Best regards,
Rookie

• This reply was modified 2 years ago by Rookie.
• This reply was modified 2 years ago by Rookie.

[Rookie]

Dear Jan,

Do you mean that I should give a case with a relatively short code? I’m sorry that I misunderstood that you think the simulation time is too long. I’m sorry that I did not succeed in writing a shorter code for you to help me, and I have really appreciated your help.I’ll try again when I have more time.

Best regards,
Rookie

[Rookie]

Dear Jan,

I believe we should focus on whether the particle simulation can run in parallel. You can even set the particle time step very small. First, you can run it in serial mode, and then switch to parallel mode to see the issue that has been bothering me.

After you run the fluid simulation without particles for the first time, in serial mode, it will directly read the fluid results fluidSolution obtained from the serial run. In parallel mode, it will directly read the fluid results fluidSolution obtained from the parallel run, so it will not simulate the fluid again.

const T physConvergeTime = 0.1 * charPhysT;                   // time until until statistics sampling in seconds
const T physStatisticsTime = 0.1 * charPhysT;                  // statistics sampling time in seconds
const T particleMaxPhysT = 0.1 * charPhysT;

Best regards,
Rookie

[Rookie]

Dear Jan,

I have successfully run it in the latest version as well. When you directly copy this code and compile it, you may encounter an error related to the particle period. In that case, replace the period code in this file with the one I sent earlier.

/olb-1.7r0/src/particles/subgrid3DLegacyFramework/boundaries/periodicBoundary3D.h

#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

Best regards,
Rookie

[Rookie]

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

[Rookie]

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 2 years, 2 months ago by Rookie.

[Rookie]

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

[Rookie]

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

[Rookie]

I followed your advice to output and read iT, and successfully implemented the continuation of the fluid simulation. I will try to write checkpoints for particles as well. Thank you for always providing me with helpful assistance. I hope the next version of OpenLB can improve this part of the functionality. Wish you all the best! I’ve placed the code below to help others solve similar problems in the future.
if (!(sLattice.load(“bstep3d.checkpoint”)))
{
for (; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, sLattice, iT, superGeometry );

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

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, converter, planeReduction, iT, superGeometry, timer );
if ( iT%( saveIter/2 )==0 && iT>0 ) {
clout << “Checkpointing the system at t=” << iT << std::endl;
sLattice.save( “bstep3d.checkpoint” );
std::ofstream checkpointFile(“bstep3d_iT.txt”);
checkpointFile << iT;
checkpointFile.close();
}
}
}
else
{
std::ifstream checkpointFile(“bstep3d_iT.txt”);
if (checkpointFile.is_open()) {
checkpointFile >> iT;
checkpointFile.close();
} else {
iT = -1;
}
for (++iT; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, sLattice, iT, superGeometry );

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

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, converter, planeReduction, iT, superGeometry, timer );
if ( iT%( saveIter/2 )==0 && iT>0 ) {
clout << “Checkpointing the system at t=” << iT << std::endl;
sLattice.save( “bstep3d.checkpoint” );
std::ofstream checkpointFile(“bstep3d_iT.txt”);
checkpointFile << iT;
checkpointFile.close();
}
}
}
timer.stop();
timer.printSummary();
}

[Rookie]

Dear Jan,

Thank you again for your response. Currently, my simulation involves first running the fluid simulation until it fully develops, and then introducing particles. I have already implemented the reverse coupling of particles to the fluid. Since the simulation time is quite long, I need to add checkpoints to prevent situations like unexpected shutdowns, ensuring that the simulation can resume from where it left off. Regarding your first question, I still haven’t found a solution. I did notice that simply uncommenting and running the checkpoint in bstep3d doesn’t successfully load and continue the simulation. I understand the importance of the if statement you emphasized, but I’ve been struggling to figure out how to incorporate it these past few days. I saw someone on the forum questioning the last time step called in the checkpoint file（https://www.openlb.net/forum/topic/checkpointing-how-to-update-it-when-loading-a-checkpoint-file/）, which raised my own doubts. Shouldn’t continuing the simulation after loading the checkpoint file mean continuing from this time step? Figuring out how to retrieve this time step seems crucial to solving the problem. I’m not sure if my modifications will achieve my goal.

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

// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp3/” );
OstreamManager clout( std::cout,”main” );
// display messages from every single mpi process
//clout.setMultiOutput(true);

UnitConverter<T,DESCRIPTOR> converter(
(T) 1./N, // physDeltaX: spacing between two lattice cells in __m__
(T) 1./(M*N), // physDeltaT: time step in __s__
(T) 1., // charPhysLength: reference length of simulation geometry
(T) 1., // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 1./100., // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1. // physDensity: physical density in __kg / m^3__
);

const int saveIter = converter.getLatticeTime( 1. );

// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
converter.write(“bstep3d”);

// === 2nd Step: Prepare Geometry ===
Vector<T,3> extend( lx0, ly0, lz0 );
Vector<T,3> origin;
IndicatorCuboid3D<T> cuboid( extend, origin );

// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 7;
#endif
CuboidGeometry3D<T> cuboidGeometry( cuboid, converter.getConversionFactorLength(), noOfCuboids );

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer );

prepareGeometry( converter, superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLattice( superGeometry );

//prepareLattice and set boundaryConditions
prepareLattice( converter, sLattice, superGeometry );

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer.start();

// Set up persistent measuring functors for result extraction
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity( sLattice, converter );
SuperEuklidNorm3D<T> normVel( velocity );

BlockReduction3D2D<T> planeReduction(
normVel,
Hyperplane3D<T>().centeredIn(cuboidGeometry.getMotherCuboid()).normalTo({0,0,1}),
600,
BlockDataSyncMode::ReduceOnly);
std::size_t iT = 0;
if (!(sLattice.load(“bstep3d.checkpoint”)))
{
for (; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, sLattice, iT, superGeometry );

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

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, converter, planeReduction, iT, superGeometry, timer );
if ( iT%( saveIter/2 )==0 && iT>0 ) {
clout << “Checkpointing the system at t=” << iT << std::endl;
sLattice.save( “bstep3d.checkpoint” );
}
}
}
else
{
for (iT = LAST_TIME_STEP_IN_THE_CHECKPOINT_FILE + 1; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, sLattice, iT, superGeometry );

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

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, converter, planeReduction, iT, superGeometry, timer );
if ( iT%( saveIter/2 )==0 && iT>0 ) {
clout << “Checkpointing the system at t=” << iT << std::endl;
sLattice.save( “bstep3d.checkpoint” );
}
}
}
timer.stop();
timer.printSummary();
}

Regarding the second issue, what I want to address is the saving and loading of SuperParticleSystem3D<T, PARTICLE> supParticleSystem(superGeometry);, because I’ve seen someone trying to implement the serialization of superGeometry（https://www.openlb.net/forum/topic/check-point-using-gpu/）, which is also achievable.I’ll be able to calculate the coupling of fluid and particles simultaneously when I stop the program and restart the computation.

Best regards,
Rookie

[Author]

Posts