[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 6 hours, 19 minutes 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

[Rookie]

I would like to express my gratitude to you and your team once again. Thank you for your previous assistance. OpenLB is indeed a fantastic software. I initially tried Palabos, but the active forum and better alignment with my research led me to choose OpenLB. I apologize for taking up your time with some simple questions (which I couldn’t figure out even after some thought).

[Rookie]

I really don’t have much knowledge about parallel processing, and I apologize for any confusion I may have caused you. The screenshot I provided earlier was from running the original channel3d program without parallel compilation but using mpirun to execute it. What I want to convey is that the large file size in the 0th step, which cannot be opened in ParaView, is due to repeated information writing, while other files can be opened normally.

I have incorporated subgrid particles into the channel3d file, and I have achieved some of the functionalities I wanted. The periodicity is implemented using this document: olb-1.6r0/src/particles/subgrid3DLegacyFramework/boundaries/periodicBoundary3D.h. I made some modifications to the original files, but these should not be the cause of the repeated writing issue. The repeated writing issue occurred only after I added some additional output. Could it be due to limitations in my computer configuration?

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);
SuperLatticePhysStrainRateFD3D<T, DESCRIPTOR> strainrate(superGeometry, sLattice, materialslist, converter);
SuperLatticePhysDissipationFD3D<T, DESCRIPTOR> dissipation(superGeometry, sLattice, materialslist, converter);
SuperLatticePhysEffectiveDissipationFD3D<T, DESCRIPTOR> effdissipation(superGeometry, sLattice, materialslist, converter, [&](Cell<T,DESCRIPTOR>& cell) -> double {
return BulkDynamics::CollisionO().computeEffectiveOmega(cell, bulkDynamicsParams);
});
SuperLatticePhysVorticityFD3D<T,DESCRIPTOR> vorticity(superGeometry, sLattice, materialslist, converter);
SuperLatticePhysStressFD3D<T,DESCRIPTOR> stress(superGeometry, sLattice, materialslist, converter);
SuperIsotropicHomogeneousTKE3D<T,DESCRIPTOR> TKE(sLattice, converter);
SuperLatticePhysEnstrophyFD3D<T,DESCRIPTOR> enstrophy(superGeometry, sLattice, materialslist, converter);
vtmWriter.addFunctor(geometry);
vtmWriter.addFunctor(velocity);
vtmWriter.addFunctor(pressure);
vtmWriter.addFunctor(sAveragedVel);
vtmWriter.addFunctor(sAveragedPre);
vtmWriter.addFunctor(sAveragedVelcc);
vtmWriter.addFunctor(wss);
vtmWriter.addFunctor(sAveragedWss);
vtmWriter.addFunctor(strainrate);
vtmWriter.addFunctor(dissipation);
vtmWriter.addFunctor(effdissipation);
vtmWriter.addFunctor(vorticity);
vtmWriter.addFunctor(stress);
vtmWriter.addFunctor(TKE);
vtmWriter.addFunctor(enstrophy);
vtmWriter.addFunctor(sAveragedVor);

vtmWriterStartTime.addFunctor(velocity);
vtmWriterStartTime.addFunctor(pressure);

[Rookie]

The screenshot I shared above is the output after I made some modifications. The following screenshot shows the output results of the original 1.5 version file, and in the zeroth step, it repeated the information.
https://postimg.cc/gallery/9TbmvYm

I added some physical quantity outputs. I believe this should not be the cause of the problem. I understand that you suggested opening MPI:
CXX :=mpic++
PARALLEL_MODE := MPI

Indeed, this allows for the correct output of VTI files. However， I need to implement periodic boundaries for particles, I can only run it with MPI turned off. Therefore, I disabled MPI, but still used “mpirun -np 8 ./channel3d” to run the program. What puzzles me is why some VTI files can be correctly output in the same program, while others show duplication.Is there any relevant setting to prevent the VTI file from repeatedly writing information after the first write?

I didn’t understand your point about completely deleting the directory. I did rename the output folder, and when outputting, two folders appear—one named “tmp” and the other with the name I specified. Is there an issue with this, and should I make changes here?

[Rookie]

Dear Adrian,

I found the case I ran before in version 1.5. At that time, I didn’t set MPI and I noticed that only the VTI file output in the 0th step was repeated 8 times, while others were normal. After modifying the files, the issue of duplicated output occurred not only in the 0th step. You can identify which ones are duplicated by observing the file sizes, and files exceeding 100MB are all duplicated. I still don’t know how to solve this problem. Then I found that this line of code is related to the setting of duplicated output:

// display messages from every single mpi process
// clout.setMultiOutput(false);

However, if this line of code is not enabled, doesn’t it mean that the output occurs only once in total for multithreading? I conducted tests on another computer, and the same issue occurred.

https://i.postimg.cc/F1t0Pq5p/repeat.png

Best regards,
Rookie

[Rookie]

Dear Adrian,

I didn’t activate MPI following the User Guide, because in a previous instance, I downloaded OpenLB-1.6 and was able to open VTK files without setting anything. Even without activating MPI, I observed 8 outputs in the terminal. After activating MPI, the terminal output occurred only once, and VTK files could be opened successfully. However, I would like to know how to ensure the correctness of my VTK files without activating MPI. Is it possible to achieve this?

Furthermore, I’ve observed that at times, among all the output files, one or several files can be opened successfully. In contrast, the files that cannot be opened exhibit duplicated outputs, and the repetition in these cases ranges from eight times to less than eight times. This also results in the files becoming very large

Best regards,
Rookie

[Rookie]

Dear Adrian

It’s not terminated, it’s when I turn off parallelism that I get these problems, it outputs the image from step 0 and the simulation doesn’t move, it doesn’t break. Normally gnuplot itself doesn’t have problems due to serial mode. I guess it could be because of this gnuplot problem or the simulation slows down when parallelism is turned off. It could be that gnuplot is not suitable for non-parallel use in this code?

Best Regards,
Rookie

[Rookie]

Dear Adrian

Yes, it outputs the image from step 0, and also as you said, this prevents the simulation from continuing. I just want to experiment more with what the code does, so I’d better remove that part of the code then.

Best Regards,
Rookie

[Rookie]

Dear Max

I’m having a similar problem, after I add these lines to my code：
vtmWriterStartTime.write(iT);
SuperEuklidNorm3D<T> normVel(velocity);
BlockReduction3D2D<T> planeReduction(normVel, {0, -1, 0}, 600, BlockDataSyncMode::ReduceOnly);
// write output as JPEG
heatmap::write(planeReduction, iT);

// write output as JPEG and changing properties
heatmap::plotParam<T> jpeg_Param;
jpeg_Param.name = “velocity”;
jpeg_Param.contourlevel = 5;
jpeg_Param.colour = “blackbody”;
jpeg_Param.zoomOrigin = {0., 0.};
jpeg_Param.zoomExtend = {1., 1.};
heatmap::write(planeReduction, iT, jpeg_Param);

if I run it in parallel it runs fine, but if I don’t I get this error:
“./channel/imageData/data/_EuklidNormphysVelocityiT0000000.p” line 21: Matrix does not represent a grid
“./channel/imageData/data/_velocityiT0000000.p” line 27: Natrix does not represent a grid

But the picture for this step is actually output.Is it that the code can only run in parallel?

Best Regards,
Rookie

[Rookie]

Dear Jan,

Thank you for the suggestion, I have downloaded this paper before, I will read this paper and other papers published by Henn, the problem I am currently experiencing is that the fluid and particles cannot be simulated at the same time, the original code simulates the fluid first and then the particles. After modifying it to reach the time step where I add the particles it reports an error. I would love to attend spring school to solve my problem, but please forgive me for not being able to attend.

terminate called after throwing an instance of ‘std::domain_error’
what(): read only access to data which is not available locally
[banana:883562] *** Process received signal ***
[banana:883562] Signal: Aborted (6)
[banana:883562] Signal code: (-6)

————————————————————————–
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
————————————————————————–
————————————————————————–
mpirun noticed that process rank 1 with PID 0 on node banana exited on signal 6 (Aborted).
————————————————————————–

Best regards,
Rookie

[Rookie]

Dear Jan,

I’m sorry to bother you again, can you point me to the literature on these two couplings, I’m trying to understand these based on the literature and the code.

simulateWithTwoWayCoupling_Mathias
simulateWithTwoWayCoupling_Davide

Best regards,
Rookie

[Author]

Posts