Particle periodicity running with MPI.

[Rookie]

Dear OpenLB developers,

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

#ifndef PERIODICBOUNDARY3D_H_
#define PERIODICBOUNDARY3D_H_

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

namespace olb {

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

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

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

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

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

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

}

#endif

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

https://postimg.cc/F7vmCHDM

Best regards,
Rookie

[jan]

Dear Rookie,

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

Here a few thoughts:

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

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

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

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

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

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

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

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

Best regards,
Jan

• This reply was modified 4 months, 2 weeks ago by jan.

[jan]

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

[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

[mathias]

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

[jan]

Dear Rookie,

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

How do you incorporate the periodic boundary? Like this?

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

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

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

Best regards,
Jan

• This reply was modified 4 months, 1 week ago by jan.

[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

[jan]

Dear Rookie,

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

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

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

Is the fluid also periodic?

Best regards,
Jan

• This reply was modified 4 months ago by jan.

[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 4 months ago by Rookie.

[jan]

Dear Rookie,

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

Best regards,
Jan

[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

[jan]

Dear Rookie,

It’s generally recommended to use the latest stable version of the software you’re working with, as it contains bug fixes and other improvements. In this case, I would like to suggest that you upgrade to the latest version of OpenLB.

In addition, it would be extremely helpful if you could provide a minimal, self-contained example that reproduces the bug you’re experiencing. This will make it easier to understand the specific circumstances leading to the problem without having to navigate through a larger codebase.

If you can provide a minimal example and confirm that you’re using the latest version of OpenLB, I’ll be better equipped to help you diagnose and fix the problem.

Best regards,
Jan

[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

[jan]

Dear Rookie,

thank you for that. With that the case indeed compiles. However, I still need a minimal example (as mentioned in a comment above) in order to help you. The provided case takes way too much time to solely compute the fluid solution. It should be possible to create a minimal example that doesn’t need any fluid calculation. Perhaps that even helps in finding the underlying problem.

Best regards,
Jan

[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

[Author]

Posts