Particle periodicity running with MPI.

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

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 3 days, 5 hours ago by 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