settlingcube3d

[yfluids]

Hello Everyone

I want to simulate sedimentation of cuboids (solid volume fraction 0.2) in fully periodic domain. I have difficulties in implementing periodic boundary conditions for particles and applying force to the fluids to ensure a constant volume flow rate. Any advice or examples would be appreciated. Thank you!

[Christoph]

Hello yFluids,
this is an adapted version of an older app. You can remove everything regarding limestones and spheres, if you only want to simulate cubes. You should also adapt boundary conditions and inflow conditions according to your needs. Note that this example can only be run with mpi enabled.

/*  Lattice Boltzmann sample, written in C++, using the OpenLB
 *  library
 *
 *  Copyright (C) 2006-2021 Nicolas Hafen, Mathias J. Krause
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  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.
 */

/*
 * TODO:
 * - Test setup with different configurations
 * - Update formatting of output as needed
 * - Use enumeration instead of preprocessor directives
 * - Improve description below
 */

/*
 * This example is based on 10.1016/j.cpc.2024.109321.
 * The limestone shapes are from the online particle database PARROT
 * (https://parrot.tu-freiberg.de/).
 * The number of triangles has been reduced.
 */

#include "olb.h"

#include <cmath>
#include <cstdint>
#include <filesystem>
#include <iostream>
#include <memory>
#include <vector>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::particles;
using namespace olb::particles::dynamics;
using namespace olb::particles::contact;
using namespace olb::particles::communication;
using namespace olb::particles::access;
using namespace olb::util;

#define WriteVTK
#define WithContact

typedef double T;

typedef enum {
    LIMESTONES,
    CUBES,
    SPHERES
} ParticleType;

const ParticleType particleType = CUBES;
T  wantedParticleVolumeFraction = 0.15;

// Define lattice type
typedef D3Q19<POROSITY, VELOCITY_NUMERATOR, VELOCITY_DENOMINATOR,
#ifdef WithContact
              CONTACT_DETECTION,
#endif
              FORCE>
    DESCRIPTOR;

// Define particle type
typedef PARTICLE_DESCRIPTOR<
    DESCRIPTOR::d, GENERAL_EXTENDABLE<DESCRIPTOR::d>,
    MOBILITY_VERLET<DESCRIPTOR::d>, SURFACE_RESOLVED_PARALLEL<DESCRIPTOR::d>,
    FORCING_RESOLVED<DESCRIPTOR::d>, PHYSPROPERTIES_RESOLVED<DESCRIPTOR::d>,
#ifdef WithContact
    MECHPROPERTIES_COLLISION<DESCRIPTOR::d>,
    NUMERICPROPERTIES_RESOLVED_CONTACT<DESCRIPTOR::d>,
#endif
    PARALLELIZATION_RESOLVED>
    PARTICLETYPE;

// Define particle-particle contact type
typedef ParticleContactArbitraryFromOverlapVolume<T, PARTICLETYPE::d, true>
    PARTICLECONTACTTYPE;

// Define particle-wall contact type
typedef WallContactArbitraryFromOverlapVolume<T, PARTICLETYPE::d, true>
    WALLCONTACTTYPE;

constexpr T nonDimensionalTime = 200.0;
constexpr T gravit             = 9.81;

T maxPhysT;

// Discretization Settings
int         res = 9;
constexpr T tau = 0.6;

// Time Settings
int iTpurge;
int iTwrite;

// Fluid Settings
constexpr T fluidDensity = 1000;
T           dynamicViscosity;
T           kinematicViscosity;

// Particle Settings
unsigned int particleNumber = 0;
T            ArchimedesNumber(1000.);
T            densityRatio(3.3);
constexpr T  radius             = 1.5e-3;
constexpr T  diameter           = 2 * radius;
T            equivalentDiameter = diameter;
T            particleDensity;
T            particleVolumeFraction;

// Domain Settings
constexpr T        extentToDiameter = 15.0;
Vector<T, 3>       extent(extentToDiameter* diameter);
const Vector<T, 3> origin(T {0});

// Characteristic Quantities
constexpr T charPhysLength = diameter;

// Contact Settings
constexpr T coefficientOfRestitution   = 0.926;
constexpr T coefficientKineticFriction = 0.16;
constexpr T coefficientStaticFriction  = T {2} * coefficientKineticFriction;
constexpr T staticKineticTransitionVelocity         = 0.001;
constexpr T YoungsModulusParticle                   = 5.0e3;
constexpr T PoissonRationParticle                   = 0.245;
T           particleEnlargementForContact           = 0;
constexpr unsigned contactBoxResolutionPerDirection = 8;

// STLs with scaling dimensions to almost reach equivalent volume
std::vector<std::pair<std::string, T>> limestoneStlFiles = {
    std::make_pair("./limestone/312_reduced.stl", 2.518e-4),
    std::make_pair("./limestone/529_reduced.stl", 2.625e-4),
    std::make_pair("./limestone/1076_reduced.stl", 1.012e-4),
    std::make_pair("./limestone/1270_reduced.stl", 1.214e-4),
    std::make_pair("./limestone/1810_reduced.stl", 0.862e-4)
  };

#ifdef PARALLEL_MODE_MPI
MPI_Comm
    averageParticleVelocityComm; /// Communicator for calculation of average particle velocity
MPI_Comm
    numberParticleCellsComm; /// Communicator for calculation of current number of particle cells
#endif // PARALLEL_MODE_MPI

std::string getParticleIdentifier(const std::size_t& pID)
{
  return std::to_string(pID);
}

T evalTerminalVelocitySingleParticleStokes()
{
  return (2. / 9.) * (particleDensity - fluidDensity) * gravit * radius *
         radius / dynamicViscosity;
};

T evalAverageSettlingVelocity(XParticleSystem<T, PARTICLETYPE>& xParticleSystem)
{
  T averageSettlingVelocity {0};

  communication::forParticlesInSuperParticleSystem<
      T, PARTICLETYPE, conditions::valid_particle_centres>(
      xParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        averageSettlingVelocity += getVelocity(particle)[2];
      });

#ifdef PARALLEL_MODE_MPI
  singleton::mpi().reduceAndBcast(averageSettlingVelocity, MPI_SUM,
                                  singleton::mpi().bossId(),
                                  averageParticleVelocityComm);
#endif // PARALLEL_MODE_MPI

  return averageSettlingVelocity / particleNumber;
}

void updateBodyForce(SuperLattice<T, DESCRIPTOR>&        sLattice,
                     XParticleSystem<T, PARTICLETYPE>&   xParticleSystem,
                     SuperGeometry<T, 3>&                superGeometry,
                     UnitConverter<T, DESCRIPTOR> const& converter)
{
  SuperLatticePhysExternalPorosity3D<T, DESCRIPTOR> porosity(sLattice,
                                                             converter);
  SuperAverage3D<T> avPorosityF(porosity, superGeometry, 1);
  int               input[1];
  T                 fluidVolumeFraction[1];
  avPorosityF(fluidVolumeFraction, input);
  const T volumeRatio = particleVolumeFraction / fluidVolumeFraction[0];

  // Apply equal to submerged weight of the particles to the fluid
  std::vector<T> balancingAcceleration(3, T(0.));
  balancingAcceleration[2] =
      gravit * volumeRatio * (particleDensity / fluidDensity - 1);

  SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
  SuperAverage3D<T>                         avgVel(velocity, superGeometry, 1);
  T                                         vel[3];
  avgVel(vel, input);
  // Avoid numerical drift
  balancingAcceleration[2] -=
      vel[2] * converter.getCharPhysVelocity() / converter.getCharPhysLength();

  const T conversionFactor = converter.getConversionFactorTime() *
                             converter.getConversionFactorTime() /
                             converter.getConversionFactorLength();
  balancingAcceleration[2] *= conversionFactor;

  AnalyticalConst3D<T, T> acc(balancingAcceleration);
  sLattice.defineField<descriptors::FORCE>(superGeometry, 1, acc);
}

template <typename T>
T calculateCubeEdgeLengthFromSphereRadius()
{
  // Calculate the volume of the sphere
  const T sphereVolume = (4.0 / 3.0) * M_PI * util::pow(radius, 3);

  // Calculate the edge length of the cube
  const T cubeEdgeLength = util::pow(sphereVolume, 1.0 / 3.0);

  return cubeEdgeLength;
}

// Prepare geometry
void prepareGeometry(SuperGeometry<T, 3>&                superGeometry,
                     UnitConverter<T, DESCRIPTOR> const& converter)
{
  OstreamManager clout(std::cout, "prepareGeometry");
  clout << "Prepare Geometry ..." << std::endl;

  superGeometry.rename(0, 1);

  Vector<T,3> origin = superGeometry.getStatistics().getMinPhysR( 1 );
  origin[1] += converter.getPhysDeltaX()/2.;
  origin[2] += converter.getPhysDeltaX()/2.;

  Vector<T,3> extend = superGeometry.getStatistics().getMaxPhysR( 1 );
  extend[0] = extend[0]-origin[0]-converter.getPhysDeltaX()/2.;
  extend[1] = extend[1]-origin[1]-converter.getPhysDeltaX()/2.;
  extend[2] = converter.getPhysDeltaX()*2.;

  IndicatorCuboid3D<T> bottom( extend, origin );
  extend[2] = converter.getPhysDeltaX()*2.;
  origin[2] = superGeometry.getStatistics().getMaxPhysR( 2 )[2]-converter.getPhysDeltaX();
  IndicatorCuboid3D<T> top (extend, origin);

  superGeometry.rename(1,3,bottom);
  superGeometry.rename(1,4,top);

  superGeometry.clean();
  superGeometry.innerClean();

  superGeometry.checkForErrors();
  superGeometry.getStatistics().print();

  clout << "Prepare Geometry ... 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)
{
  OstreamManager clout(std::cout, "prepareLattice");
  clout << "Prepare Lattice ..." << std::endl;

  /// Material=0 -->do nothing
  sLattice.defineDynamics<
      PorousParticleKupershtokhForcedBGKdynamics<T, DESCRIPTOR>>(superGeometry,
                                                                 1);
  sLattice.defineDynamics<BounceBack>(superGeometry, 2);

  sLattice.setParameter<descriptors::OMEGA>(
      converter.getLatticeRelaxationFrequency());

  {
    auto& communicator = sLattice.getCommunicator(stage::PostPostProcess());
    communicator
        .requestFields<POROSITY, VELOCITY_NUMERATOR, VELOCITY_DENOMINATOR>();
    communicator.requestOverlap(sLattice.getOverlap());
    communicator.exchangeRequests();
  }
  

  clout << "Prepare Lattice ... OK" << std::endl;
}

// Set Boundary Values
void setBoundaryValues(SuperLattice<T, DESCRIPTOR>&        sLattice,
                       UnitConverter<T, DESCRIPTOR> const& converter, int iT,
                       SuperGeometry<T, 3>& superGeometry)
{
  OstreamManager clout(std::cout, "setBoundaryValues");

  boundary::set<boundary::LocalVelocity>(sLattice, superGeometry, 3); //inlet
  boundary::set<boundary::LocalPressure>(sLattice, superGeometry, 4); //outlet

  if (iT == 0) {
    AnalyticalConst3D<T, T> zero(0.);
    AnalyticalConst3D<T, T> one(1.);
    sLattice.defineField<descriptors::POROSITY>(
        superGeometry.getMaterialIndicator(1), one);
    // Set initial condition
    AnalyticalConst3D<T, T>    ux(0.);
    AnalyticalConst3D<T, T>    uy(0.);
    AnalyticalConst3D<T, T>    uz(0.);
    AnalyticalComposed3D<T, T> u(ux, uy, uz);

    AnalyticalConst3D<T, T> rho(1.);

    // Initialize all values of distribution functions to their local equilibrium
    sLattice.defineRhoU(superGeometry, 1, rho, u);
    sLattice.iniEquilibrium(superGeometry, 1, rho, u);

    // Make the lattice ready for simulation
    sLattice.initialize();
  }
}

void getResults(SuperLattice<T, DESCRIPTOR>&        sLattice,
                UnitConverter<T, DESCRIPTOR> const& converter, int iT,
                SuperGeometry<T, 3>& superGeometry, Timer<double>& timer,
                XParticleSystem<T, PARTICLETYPE>& xParticleSystem)
{
  OstreamManager clout(std::cout, "getResults");

  SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
#ifdef WriteVTK
  SuperVTMwriter3D<T>                       vtkWriter("sedimentation");
  SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
  SuperLatticePhysExternalPorosity3D<T, DESCRIPTOR> externalPor(sLattice,
                                                                converter);

  vtkWriter.addFunctor(pressure);
  vtkWriter.addFunctor(velocity);
  vtkWriter.addFunctor(externalPor);

  if (iT == 0) {
    /// Writes the converter log file
    SuperLatticeCuboid3D<T, DESCRIPTOR>   cuboid(sLattice);
    SuperLatticeRank3D<T, DESCRIPTOR>     rank(sLattice);
    vtkWriter.write(cuboid);
    vtkWriter.write(rank);
    vtkWriter.createMasterFile();
  }

  if (iT % iTwrite == 0) {
    vtkWriter.write(iT);
  }
#endif // WriteVTK

  /// Writes output on the console
  if (iT % iTwrite == 0) {
    // TODO: Update formatting as needed
    clout << "Average settling velocity: "
          << evalAverageSettlingVelocity(xParticleSystem) << " in m/s"
          << std::endl;

    timer.update(iT);
    timer.printStep();
    sLattice.getStatistics().print(iT, converter.getPhysTime(iT));
  }
}

int main(int argc, char* argv[])
{
  /// === 1st Step: Initialization ===
  initialize(&argc, &argv);
  OstreamManager clout(std::cout, "main");
  #ifdef PARALLEL_MODE_MPI //Check if MPI is available, otherwise throw error
  std::vector<std::string> cmdInput;

  if (argc > 1) {
    cmdInput.assign(argv + 1, argv + argc);
    res = std::stoi(cmdInput[0]);
  }
  if (argc > 2) {
    wantedParticleVolumeFraction = std::stod(cmdInput[1]);
  }
  if (argc > 3) {
    densityRatio = std::stod(cmdInput[2]);
  }
  if (argc > 4) {
    ArchimedesNumber = std::stod(cmdInput[3]);
  }

  const std::string positionsfilename =[]{
    if constexpr (particleType==LIMESTONES){
      return std::string("limestone/particlepositions_limestone_15_0.150021");
    }
    else return std::string("particlepositions_sphere_1.5mm_15_0.30");
  }();

  const T domainVolume = extent[0] * extent[1] * extent[2];
  particleDensity      = densityRatio * fluidDensity;
  kinematicViscosity =
      util::sqrt(gravit * (densityRatio - 1) *
                 util::pow(equivalentDiameter, 3) / ArchimedesNumber);

  const Vector<T, 3> externalAcceleration = {
      .0, .0, -gravit * (1. - fluidDensity / particleDensity)};

  // Estimation maximal velocity using Stoke's law
  dynamicViscosity         = kinematicViscosity * fluidDensity;
  const T charPhysVelocity = evalTerminalVelocitySingleParticleStokes();
  UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
      (int)res,              // resolution
      (T)tau,                // latticeRelaxationTime
      (T)charPhysLength,     // charPhysLength
      (T)charPhysVelocity,   // charPhysVelocity
      (T)kinematicViscosity, // physViscosity
      (T)fluidDensity        // fluidDensity
  );
  converter.write();
  converter.print();

  if (MPI_Comm_dup(MPI_COMM_WORLD, &averageParticleVelocityComm) !=
      MPI_SUCCESS) {
    throw std::runtime_error("Unable to duplicate MPI communicator");
  }
  if (MPI_Comm_dup(MPI_COMM_WORLD, &numberParticleCellsComm) != MPI_SUCCESS) {
    throw std::runtime_error("Unable to duplicate MPI communicator");
  }

  std::size_t iT = 0;
  maxPhysT       = nonDimensionalTime * equivalentDiameter / charPhysVelocity;
  particleEnlargementForContact = converter.getPhysDeltaX() / T {5};

  /// === 2rd Step: Prepare Geometry ===
  /// Instantiation of a cuboidGeometry with weights
  IndicatorCuboid3D<T> cuboid(extent, origin);
  constexpr auto       getPeriodicity = []() {
    return Vector<bool, 3>(true, true, true);
  };

  const unsigned numberProcesses = singleton::mpi().getSize();

CuboidDecomposition3D<T> cuboidGeometry(
      cuboid, converter.getConversionFactorLength(), numberProcesses);

cuboidGeometry.setPeriodicity({ true, false, false });

  HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
  SuperGeometry<T, 3>      superGeometry(cuboidGeometry, loadBalancer, 2);
  prepareGeometry(superGeometry, converter);

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

  prepareLattice(sLattice, converter, superGeometry);

  std::vector<SolidBoundary<T, DESCRIPTOR::d>> solidBoundaries;
  const T epsilon     = T {0.5} * converter.getConversionFactorLength();
  const T halfEpsilon = T {0.5} * epsilon;

  constexpr T overlapSecurityFactor = 1.0;
  T maxCircumRadius = T {0};

  // Create smooth indicators for limestones
  std::vector<std::shared_ptr<STLreader<T>>> limestoneSTLreaders;
  std::vector<std::unique_ptr<olb::SmoothIndicatorCustom3D<T, T, true>>>
      limestoneIndicators;
  const T latticeSpacingDiscreteParticle =
    T {0.2} * converter.getConversionFactorLength();

  if constexpr (particleType==LIMESTONES){
    clout << "Initializing limestone ..." << std::endl;

    {
      unsigned i = 0;
      for (auto& STLfile : limestoneStlFiles) {
        limestoneSTLreaders.push_back(std::make_shared<STLreader<T>>(
            STLfile.first, converter.getConversionFactorLength(),
            STLfile.second));
        limestoneIndicators.push_back(
            std::make_unique<olb::SmoothIndicatorCustom3D<T, T, true>>(
                latticeSpacingDiscreteParticle, limestoneSTLreaders[i],
                olb::Vector<T, 3>(T {}), epsilon, olb::Vector<T, 3>(T {})));
        limestoneIndicators.back()->calcMofiAndMass(particleDensity);
        ++i;
      }
    }

    clout << "Initializing limestone ... OK" << std::endl;

    // Create Particle Dynamics
    // Create ParticleSystem

    for (auto& STLsurface : limestoneIndicators) {
      maxCircumRadius = util::max(maxCircumRadius, STLsurface->getCircumRadius());
    }
  }
  if constexpr (particleType==SPHERES){
    maxCircumRadius = radius + halfEpsilon;
  }

  if constexpr (access::providesContactMaterial<PARTICLETYPE>()) {
    const T detectionDistance =
        T {0.5} * util::sqrt(PARTICLETYPE::d) * converter.getPhysDeltaX();
    maxCircumRadius = maxCircumRadius - halfEpsilon +
                      util::max(halfEpsilon, detectionDistance);
  }
  maxCircumRadius *= overlapSecurityFactor;

  // ensure parallel mode is enabled

    SuperParticleSystem<T, PARTICLETYPE> xParticleSystem(superGeometry,
                                                       maxCircumRadius);
    particles::communication::checkCuboidSizes(xParticleSystem);

  // Create ParticleContact container
  ContactContainer<T, PARTICLECONTACTTYPE, WALLCONTACTTYPE> contactContainer;
  // Container containg contact properties
  ContactProperties<T, 1> contactProperties;
  contactProperties.set(
      0, 0,
      evalEffectiveYoungModulus(YoungsModulusParticle, YoungsModulusParticle,
                                PoissonRationParticle, PoissonRationParticle),
      coefficientOfRestitution, coefficientKineticFriction,
      coefficientStaticFriction, staticKineticTransitionVelocity);

  // Create and assign resolved particle dynamics
  xParticleSystem.defineDynamics<VerletParticleDynamics<T, PARTICLETYPE>>();

  //Create particle manager handling coupling, gravity and particle dynamics
  ParticleManager<T, DESCRIPTOR, PARTICLETYPE> particleManager(
      xParticleSystem, superGeometry, sLattice, converter, externalAcceleration, getPeriodicity());

  // Create Communicators

  const communication::ParticleCommunicator& particleCommunicator =
      particleManager.getParticleCommunicator();

  const std::function<T(const std::size_t&)> getCircumRadius =
      [&](const std::size_t& pID) {
        if constexpr (particleType==LIMESTONES){
          return limestoneIndicators[pID % limestoneStlFiles.size()]
              ->getCircumRadius();
        }
        if constexpr (particleType==SPHERES){
          return radius + T {0.5} * epsilon;
        }
        if constexpr (particleType==CUBES){

          return T {0.5} *
                (calculateCubeEdgeLengthFromSphereRadius<T>() * util::sqrt(3) +
                  epsilon);
        }

      };
  const std::function<T(const std::size_t&)> getParticleVolume =
      [&](const std::size_t& pID) {
        if constexpr (particleType==LIMESTONES){
          return limestoneIndicators[pID % limestoneStlFiles.size()]->getVolume();
        }
        if constexpr (particleType==SPHERES){
          return T {4} / T {3} * M_PI * radius * radius * radius;
        }
        if constexpr (particleType==CUBES){
          return util::pow(calculateCubeEdgeLengthFromSphereRadius<T>(), 3);
        }

      };

  const std::function<void(
      const particles::creators::SpawnData<T, DESCRIPTOR::d>&,
      const std::string&)>
      createParticleFromString =
          [&](const particles::creators::SpawnData<T, DESCRIPTOR::d>& data,
              const std::string&                                      pID) {
            const PhysR<T, 3>  physPosition   = data.position;
            const Vector<T, 3> angleInDegrees = data.angleInDegree;

          if constexpr (particleType==LIMESTONES){
            if (particleNumber < limestoneStlFiles.size()) {
              std::shared_ptr <STLreader<T>> limestoneIndicator =
                  std::make_shared<STLreader<T>>(
                      limestoneStlFiles[particleNumber].first,
                      converter.getConversionFactorLength(),
                      limestoneStlFiles[particleNumber].second);
              creators::addResolvedArbitraryShape3D<T, PARTICLETYPE>(
                  xParticleSystem,  physPosition,
                  latticeSpacingDiscreteParticle,
                  limestoneIndicator,
                  epsilon, particleDensity );
            }
            else {
              creators::addResolvedObject<T, PARTICLETYPE>(
                  xParticleSystem,
                  particleNumber % limestoneStlFiles.size(), physPosition,
                  particleDensity, angleInDegrees);
            }
          }
          if constexpr (particleType==SPHERES){
            creators::addResolvedSphere3D(xParticleSystem,
                                          physPosition, radius, epsilon,
                                          particleDensity);
          }
          if constexpr (particleType==CUBES){
            creators::addResolvedCuboid3D(
                xParticleSystem, physPosition,
                Vector<T, 3>(calculateCubeEdgeLengthFromSphereRadius<T>()),
                epsilon, particleDensity);
          }
            ++particleNumber;
          };

  if (positionsfilename != "") {
    clout << "Spawning particles from " << positionsfilename << " ..."
          << std::endl;
    if (std::filesystem::exists(positionsfilename)) {
      std::vector<particles::creators::SpawnData<T, DESCRIPTOR::d>>
          tmpSpawnData;
          tmpSpawnData = particles::creators::setParticles<T, 3>(
          positionsfilename, wantedParticleVolumeFraction, cuboid, domainVolume,
          getParticleVolume, createParticleFromString);

      T tmpVolume = T {0};
      for (unsigned pID = 0; pID < tmpSpawnData.size(); ++pID) {
        tmpVolume += getParticleVolume(pID);
      }
      particleVolumeFraction = tmpVolume / domainVolume;
    }
    else {
      OLB_ASSERT(false, positionsfilename + " does not exist.");
    }
  }
  else {
    OLB_ASSERT(false, "No particle positions file given.");
  }
 if constexpr (particleType==LIMESTONES){
    limestoneIndicators.clear();
  }

  clout << "Spawning particles from " << positionsfilename << " ... OK"
        << std::endl;

  maxCircumRadius = 0.;
  forParticlesInSuperParticleSystem(
      xParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
#ifdef WithContact
        particle.setField<MECHPROPERTIES, MATERIAL>(0);
        particle.setField<NUMERICPROPERTIES, ENLARGEMENT_FOR_CONTACT>(
            particleEnlargementForContact);
#endif // WithContact

        const T currCircumRadius = access::getRadius(particle);
        maxCircumRadius          = util::max(currCircumRadius, maxCircumRadius);
      });

  singleton::mpi().reduceAndBcast(maxCircumRadius, MPI_MAX,
                                  singleton::mpi().bossId(),
                                  particleCommunicator.contactTreatmentComm);

  xParticleSystem.updateOffsetFromCircumRadius(overlapSecurityFactor *
                                               maxCircumRadius);

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

  /// === 4th Step: Main Loop with Timer ===
  clout << "MaxIT: " << converter.getLatticeTime(maxPhysT) << std::endl;
  iTwrite = util::max(0.02 * converter.getLatticeTime(maxPhysT), 1);
  iTpurge = util::max(util::ceil(0.06 * converter.getLatticeTime(maxPhysT)), 1);

  Timer<double> timer(converter.getLatticeTime(maxPhysT),
                      superGeometry.getStatistics().getNvoxel());
  timer.start();
  for (iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) {
#ifndef WithContact
    // Execute particle manager
    particleManager.execute<
        couple_lattice_to_parallel_particles<T, DESCRIPTOR, PARTICLETYPE>,
        communicate_parallel_surface_force<T, PARTICLETYPE>,
        apply_gravity<T, PARTICLETYPE>,
        process_dynamics_parallel<T, PARTICLETYPE>,
        update_particle_core_distribution<T, PARTICLETYPE>>();

    particles::dynamics::coupleResolvedParticlesToLattice<
        T, DESCRIPTOR, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE>(
        xParticleSystem, contactContainer, superGeometry, sLattice, converter,
        solidBoundaries, getPeriodicity);

#else  // WithContact
    // Couple lattice to particles
    couple_lattice_to_parallel_particles<T, DESCRIPTOR, PARTICLETYPE>::execute(
        xParticleSystem, superGeometry, sLattice, converter, getPeriodicity());

    communicate_parallel_surface_force<T, PARTICLETYPE>::execute(
        xParticleSystem, particleCommunicator);

    // Apply contacts
    particles::contact::processContacts<T, PARTICLETYPE, PARTICLECONTACTTYPE,
                                        WALLCONTACTTYPE,
                                        ContactProperties<T, 1>>(
        xParticleSystem, solidBoundaries, contactContainer, contactProperties,
        superGeometry,
        particleCommunicator.contactTreatmentComm,
        contactBoxResolutionPerDirection, T {4. / (3 * util::sqrt(M_PI))},
        getPeriodicity);

    // Apply gravity
    communication::forParticlesInSuperParticleSystem<
        T, PARTICLETYPE,
        conditions::valid_particle_centres //only consider center for resolved
        >(xParticleSystem,
          [&](Particle<T, PARTICLETYPE>&       particle,
              ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
            apply_gravity<T, PARTICLETYPE>::execute(xParticleSystem, particle,
                                                    externalAcceleration,
                                                    converter.getPhysDeltaT());
          });

    // Process particles (Verlet algorithm)
    communication::forParticlesInSuperParticleSystem<
        T, PARTICLETYPE,
        conditions::valid_particle_centres //only consider center for resolved
        >(xParticleSystem,
          [&](Particle<T, PARTICLETYPE>&       particle,
              ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
            particle.process(converter.getPhysDeltaT());
          });

    communicatePostContactTreatmentContacts(
        contactContainer, xParticleSystem, converter.getPhysDeltaX(),
        particleCommunicator.particleContactDetectionComm,
        particleCommunicator.wallContactDetectionComm,
        getPeriodicity());

    update_particle_core_distribution<T, PARTICLETYPE>::execute(
        xParticleSystem, converter.getPhysDeltaX(), particleCommunicator,
        getPeriodicity());

    // Couple particles to lattice
    particles::dynamics::coupleResolvedParticlesToLattice<
        T, DESCRIPTOR, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE>(
        xParticleSystem, contactContainer, superGeometry, sLattice, converter,
        solidBoundaries, getPeriodicity);

    // Communicate found contacts
    communicateContactsParallel(
        contactContainer, xParticleSystem, converter.getPhysDeltaX(),
        particleCommunicator.particleContactDetectionComm,
        particleCommunicator.wallContactDetectionComm,
        getPeriodicity());

    if constexpr (isPeriodic(getPeriodicity())) {
      accountForPeriodicParticleBoundary(xParticleSystem, contactContainer,
                                         superGeometry, getPeriodicity);
    }
#endif // WithContact

    if (iT % iTpurge == 0) {
      purgeInvalidParticles<T, PARTICLETYPE>(xParticleSystem);
    }

    // Get Results
    getResults(sLattice, converter, iT, superGeometry, timer, xParticleSystem);

    updateBodyForce(sLattice, xParticleSystem, superGeometry, converter);

    // Collide and stream
    sLattice.collideAndStream();
  }

  timer.update(iT);
  timer.stop();
  timer.printSummary();
  #else
    throw std::runtime_error(
      "This example is not designed to run in serial mode.");
  #endif // PARALLEL_MODE_MPI
}

If you have more questions, feel free to ask for more guidance.

Kind regards,
Christoph

[yfluids]

Hi Christoph,

Thanks very much. When I ran the code using mpirun, it did not work. Could you please tell me What is the command to run the simulation? I noticed that the particle creator function in this code is different from the function in particleCreatorFunctions3D.h. The orientation of the particle is not specified. How is the orientation described in the new version code?

Thank you.

[Christoph]

Here is a modified and drastically simplified version of the example. This does compile in the current release directory.
Please modify to match your required volume fraction and add inlet and outlet.
If you have any more questions please ask again.
With kind regards,
Christoph
´
/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2006-2021 Nicolas Hafen, 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.
*/

/* settlingCube3d.cpp:
* The case examines the settling of a cubical silica particle
* under the influence of gravity.
* The object is surrounded by water in a rectangular domain
* limited by no-slip boundary conditions.
* For the calculation of forces an DNS approach is chosen
* which also leads to a back-coupling of the particle on the fluid,
* inducing a flow.
*
* The simulation is based on the homogenised lattice Boltzmann approach
* (HLBM) introduced in “Particle flow simulations with homogenised
* lattice Boltzmann methods” by Krause et al.
* and extended in “Towards the simulation of arbitrarily shaped 3D particles
* using a homogenised lattice Boltzmann method” by Trunk et al.
* for the simulation of 3D particles.
*
* This example demonstrates the usage of HLBM in the OpenLB framework.
* To improve parallel performance, the particle decomposition scheme
* described in 10.48550/arXiv.2312.14172 is used.
*/

#include “olb3D.h”
#include “olb3D.hh” // use generic version only!
#include <random>
#include <vector>
#define random

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::particles;
using namespace olb::particles::dynamics;
using namespace olb::particles::contact;
using namespace olb::particles::communication;
using namespace olb::particles::access;
using namespace olb::util;
typedef double T;
// Define lattice type
typedef PorousParticleD3Q19Descriptor DESCRIPTOR;

// Define particleType
#ifdef PARALLEL_MODE_MPI
// Particle decomposition improves parallel performance
typedef ResolvedDecomposedParticle3D PARTICLETYPE;
#else
typedef ResolvedParticle3D PARTICLETYPE;
#endif
// Define particle-particle contact type
typedef ParticleContactArbitraryFromOverlapVolume<T, PARTICLETYPE::d, true>
PARTICLECONTACTTYPE;
// Define particle-wall contact type
typedef WallContactArbitraryFromOverlapVolume<T, PARTICLETYPE::d, true>
WALLCONTACTTYPE;

#define WriteVTK

// Discretization Settings
int res = 30;
T const charLatticeVelocity = 0.01;

// Time Settings
T const maxPhysT = 0.8; // max. simulation time in s
T const iTwrite = 0.002; // write out intervall in s

// Domain Settings
T const lengthX = 0.05;
T const lengthY = 0.05;
T const lengthZ = 0.05;

// Fluid Settings
T const physDensity = 1000;
T const physViscosity = 1E-5;

//Particle Settings
int numParticles = 11;
T centerX = lengthX*.5;
T centerY = lengthY*.5;
T centerZ = lengthZ*.5;
T const cubeDensity = 2500;
T const cubeEdgeLength = 0.0025;
T const smallCubeEdgeLength = 0.0015;
Vector<T,3> cubeCenter = {centerX,centerY,centerZ};
Vector<T,3> cubeOrientation = {0.,15.,0.};
Vector<T,3> cubeVelocity = {0.,0.,0.};
Vector<T,3> externalAcceleration = {.0, .0, -T(9.81) * (T(1) – physDensity / cubeDensity)};

// Characteristic Quantities
T const charPhysLength = lengthX;
T const charPhysVelocity = 0.15; // Assumed maximal velocity

// Prepare geometry
void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry )
{
OstreamManager clout(std::cout, “prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 1);

superGeometry.clean();
superGeometry.innerClean();

superGeometry.checkForErrors();
superGeometry.getStatistics().print();
clout << “Prepare Geometry … OK” << std::endl;
return;
}

// Set up the geometry of the simulation
void prepareLattice(
SuperLattice<T, DESCRIPTOR>& sLattice, UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “prepareLattice”);
clout << “Prepare Lattice …” << std::endl;
clout << “setting Velocity Boundaries …” << std::endl;

sLattice.defineDynamics<PorousParticleBGKdynamics>(superGeometry, 1); //fluid

sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());

{
auto& communicator = sLattice.getCommunicator(stage::PostPostProcess());
communicator.requestFields<POROSITY,VELOCITY_NUMERATOR,VELOCITY_DENOMINATOR>();
communicator.requestOverlap(sLattice.getOverlap());
communicator.exchangeRequests();
}

clout << “Prepare Lattice … OK” << std::endl;
}

//Set Boundary Values
void setBoundaryValues(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “setBoundaryValues”);

if (iT == 0) {
AnalyticalConst3D<T, T> zero(0.);
AnalyticalConst3D<T, T> one(1.);
sLattice.defineField<descriptors::POROSITY>(superGeometry.getMaterialIndicator({0,1,2}), one);

// Set initial condition
AnalyticalConst3D<T, T> ux(0.);
AnalyticalConst3D<T, T> uy(0.);
AnalyticalConst3D<T, T> uz(0.);
AnalyticalConst3D<T, T> rho(1.);
AnalyticalComposed3D<T, T> u(ux, uy, uz);

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU(superGeometry, 1, rho, u);
sLattice.iniEquilibrium(superGeometry, 1, rho, u);

// 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
void getResults(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry, Timer<T>& timer,
XParticleSystem<T, PARTICLETYPE>& xParticleSystem )
{
OstreamManager clout(std::cout, “getResults”);

#ifdef WriteVTK
SuperVTMwriter3D<T> vtkWriter(“sedimentation”);
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
SuperLatticePhysExternalPorosity3D<T, DESCRIPTOR> externalPor(sLattice, converter);
vtkWriter.addFunctor(velocity);
vtkWriter.addFunctor(pressure);
vtkWriter.addFunctor(externalPor);

if (iT == 0) {
/// Writes the converter log file
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}

if (iT % converter.getLatticeTime(iTwrite) == 0) {
vtkWriter.write(iT);
}
#endif

/// Writes output on the console
if (iT % converter.getLatticeTime(iTwrite) == 0) {
timer.update(iT);
timer.printStep();
sLattice.getStatistics().print(iT, converter.getPhysTime(iT));
#ifdef PARALLEL_MODE_MPI
communication::forParticlesInSuperParticleSystem<
T, PARTICLETYPE, conditions::valid_particle_centres>(
xParticleSystem,
[&](Particle<T, PARTICLETYPE>& particle,
ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
io::printResolvedParticleInfo(particle);
});
#else
for (std::size_t iP=0; iP<xParticleSystem.size(); ++iP) {
auto particle = xParticleSystem.get(iP);
io::printResolvedParticleInfo(particle);
}
#endif
}
}

int main(int argc, char* argv[])
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout, “main”);

UnitConverterFromResolutionAndLatticeVelocity<T,DESCRIPTOR> converter(
(int) res, //resolution
( T ) charLatticeVelocity, //charLatticeVelocity
( T ) charPhysLength, //charPhysLength
( T ) charPhysVelocity, //charPhysVelocity
( T ) physViscosity, //physViscosity
( T ) physDensity //physDensity
);
converter.print();

/// === 2rd Step: Prepare Geometry ===
/// Instantiation of a cuboidGeometry with weights
Vector<T,3> origin( 0. );
Vector<T,3> extend( lengthX, lengthY, lengthZ );
IndicatorCuboid3D<T> cuboid(extend, origin);

// std::string fName(“sedimentation.xml”);
// XMLreader config(fName);

#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());
#else
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), 7);
#endif
cuboidGeometry.print();

cuboidGeometry.setPeriodicity( true, true, true );

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T,3> superGeometry(cuboidGeometry, loadBalancer, 2);

prepareGeometry(converter,superGeometry);

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

// Prepare lattice
prepareLattice(sLattice, converter, superGeometry);

// Create ParticleSystem
#ifdef PARALLEL_MODE_MPI
SuperParticleSystem<T,PARTICLETYPE> particleSystem(superGeometry);
#else
ParticleSystem<T,PARTICLETYPE> particleSystem;
#endif

//Create particle manager handling coupling, gravity and particle dynamics
ParticleManager<T,DESCRIPTOR,PARTICLETYPE> particleManager(
particleSystem, superGeometry, sLattice, converter, externalAcceleration);

// Create and assign resolved particle dynamics
particleSystem.defineDynamics<
VerletParticleDynamics<T,PARTICLETYPE>>();

// Calculate particle quantities
T epsilon = 0.5*converter.getPhysDeltaX();
Vector<T,3> cubeExtend( cubeEdgeLength );
Vector<T,3> smallCubeExtend( smallCubeEdgeLength );

T min= lengthX*T(0.2);
T max= lengthX*T(0.8);
std::vector<T> positionsX;
std::vector<T> positionsY;
std::vector<T> positionsZ;

for(int i=1; i<numParticles; i++){
positionsX.push_back(min+(max-min)*i/numParticles);
positionsY.push_back(min+(max-min)*i/numParticles);
positionsZ.push_back(min+(max-min)*i/numParticles);
}

for( int i=1; i<numParticles ; i++ ){
for( int j=1; j<numParticles; j++ ){
for( int k=1; k<numParticles; k++ ){
if( positionsX[i] > min && positionsX[i] < max &&
positionsY[j] > min && positionsY[j] < max &&
positionsZ[k] > min && positionsZ[k] < max ){
cubeCenter = { positionsX[i] , positionsY[j] ,positionsZ[k] };
creators::addResolvedCuboid3D( particleSystem, cubeCenter,
smallCubeExtend, epsilon, cubeDensity, cubeOrientation );
}
}
}
}

// Check ParticleSystem
particleSystem.checkForErrors();

/// === 4th Step: Main Loop with Timer ===
Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel());
timer.start();

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

clout << “MaxIT: ” << converter.getLatticeTime(maxPhysT) << std::endl;

for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT)+10; ++iT) {

// Execute particle manager
particleManager.execute<
couple_lattice_to_particles<T,DESCRIPTOR,PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
communicate_surface_force<T,PARTICLETYPE>,
#endif
apply_gravity<T,PARTICLETYPE>,
process_dynamics<T,PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
update_particle_core_distribution<T, PARTICLETYPE>,
#endif
couple_particles_to_lattice<T,DESCRIPTOR,PARTICLETYPE>
>();

// Get Results
getResults(sLattice, converter, iT, superGeometry, timer, particleSystem );

// Collide and stream
sLattice.collideAndStream();
}

timer.stop();
timer.printSummary();
}
`

[2366039968@qq.com]

Hello everyone, I modified the settlingCube3d code to simulate an oblique collision between a particle and a wall. I want to run this modified code on the GPU. My environment has already been configured properly, and the original code can run on the GPU without any issues. However, when I run my modified version in the GPU environment, the GPU memory usage stays at zero. I’m not sure if this is caused by incompatibility in my code.

Below is the code after my modifications.

/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2006-2021 Nicolas Hafen, 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.
*/

/* settlingCube3d.cpp:
* The case examines the settling of a cubical silica particle
* under the influence of gravity.
* The object is surrounded by water in a rectangular domain
* limited by no-slip boundary conditions.
* For the calculation of forces an DNS approach is chosen
* which also leads to a back-coupling of the particle on the fluid,
* inducing a flow.
*
* The simulation is based on the homogenised lattice Boltzmann approach
* (HLBM) introduced in “Particle flow simulations with homogenised
* lattice Boltzmann methods” by Krause et al.
* and extended in “Towards the simulation of arbitrarily shaped 3D particles
* using a homogenised lattice Boltzmann method” by Trunk et al.
* for the simulation of 3D particles.
*
* This example demonstrates the usage of HLBM in the OpenLB framework.
* To improve parallel performance, the particle decomposition scheme
* described in 10.48550/arXiv.2312.14172 is used.
*/

#include “olb3D.h”
#include “olb3D.hh” // use generic version only!

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace olb::particles;
using namespace olb::particles::dynamics;
using namespace std;

using T = FLOATING_POINT_TYPE;

// Define lattice type
typedef PorousParticleD3Q19Descriptor DESCRIPTOR;

// Define particleType
#ifdef PARALLEL_MODE_MPI
// Particle decomposition improves parallel performance
typedef ResolvedDecomposedParticle3D PARTICLETYPE;
#else
typedef ResolvedParticle3D PARTICLETYPE;
#endif

#define WriteVTK
#define WriteGnuPlot

//std::string gnuplotFilename = “gnuplot.dat”;
// Discretization Settings
int res = 800;
T const charLatticeVelocity =0.03;

// Time Settings
T const maxPhysT = 0.15; // max. simulation time in s
T const iTwrite = 0.001; // write out intervall in s

// Domain Settings
T const lengthX = 0.1016;
T const lengthY = 0.1016;
T const lengthZ = 0.1016;

// Fluid Settings
T const physDensity = 1000;
T const physViscosity = 1E-6;
T const sphered = 0.00635/2;
//Particle Settings
T centerX = lengthX*.5;
T centerY =0.0079268;
T centerZ = lengthZ*.5;
T const sphereDensity = 4000;

Vector<T,3> sphereCenter = {centerX,centerY,centerZ};
// External acceleration: gravity
Vector<T,3> externalAcceleration = {.0, .0, -9.81};

// Particle settings: ensure initial velocity is zero
Vector<T,3> sphereVelocity = {0., -0.07884, 0.07884};

// Characteristic Quantities./s
T const charPhysLength = lengthX;
T const charPhysVelocity = 0.2; // Assumed maximal velocity

// Prepare geometry
void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 2);
superGeometry.rename(2, 1, {1, 1, 1});

superGeometry.clean();
superGeometry.innerClean();

superGeometry.checkForErrors();
superGeometry.getStatistics().print();
clout << “Prepare Geometry … OK” << std::endl;
return;
}

// Set up the geometry of the simulation
void prepareLattice(
SuperLatticeGpu<T, DESCRIPTOR> sLattice(superGeometry);
, UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “prepareLattice”); //output information control
clout << “Prepare Lattice …” << std::endl; //out prepare information
clout << “setting Velocity Boundaries …” << std::endl; //velocity boundaries

/// Material=0 –>do nothing
sLattice.defineDynamics<PorousParticleBGKdynamics>(superGeometry, 1); //dynamics module
setBounceBackBoundary(sLattice, superGeometry, 2); // bounce boundary

sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());//relax OMEGA

{
auto& communicator = sLattice.getCommunicator(stage::PostPostProcess());
communicator.requestFields<POROSITY,VELOCITY_NUMERATOR,VELOCITY_DENOMINATOR>();//kongxilv sudufenzi,fenmu
communicator.requestOverlap(sLattice.getOverlap());//overlap area
communicator.exchangeRequests();
}

clout << “Prepare Lattice … OK” << std::endl; //papare ok
}

//Set Boundary Values
void setBoundaryValues(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “setBoundaryValues”);

if (iT == 0) {
AnalyticalConst3D<T, T> zero(0.);
AnalyticalConst3D<T, T> one(1.);
sLattice.defineField<descriptors::POROSITY>(superGeometry.getMaterialIndicator({0,1,2}), one);
// Set initial condition
AnalyticalConst3D<T, T> ux(0.);
AnalyticalConst3D<T, T> uy(0.);
AnalyticalConst3D<T, T> uz(0.);
AnalyticalConst3D<T, T> rho(1.);
AnalyticalComposed3D<T, T> u(ux, uy, uz);

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU(superGeometry, 1, rho, u);
sLattice.iniEquilibrium(superGeometry, 1, rho, u);

// Make the lattice ready for simulation
sLattice.initialize();
}
}

/// Computes the pressure drop between the voxels before and after the cylinder
void getResults(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry, Timer<T>& timer,
XParticleSystem<T, PARTICLETYPE>& xParticleSystem )
{
OstreamManager clout(std::cout, “getResults”); //information control

#ifdef WriteVTK
SuperVTMwriter3D<T> vtkWriter(“sedimentation”);
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
SuperLatticePhysExternalPorosity3D<T, DESCRIPTOR> externalPor(sLattice, converter);
vtkWriter.addFunctor(velocity);
vtkWriter.addFunctor(pressure);
vtkWriter.addFunctor(externalPor);

if (iT == 0) {
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}

if (iT % converter.getLatticeTime(iTwrite) == 0) {
vtkWriter.write(iT);
}
#endif

#ifdef PARALLEL_MODE_MPI
communication::forParticlesInSuperParticleSystem<
T, PARTICLETYPE, conditions::valid_particle_centres>(
xParticleSystem,
[&](Particle<T, PARTICLETYPE>& particle,
ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
io::printResolvedParticleInfo(particle);
});
#else
for (std::size_t iP=0; iP<xParticleSystem.size(); ++iP) {
auto particle = xParticleSystem.get(iP);
// io::printResolvedParticleInfo(particle);
}
#endif

}

int main(int argc, char* argv[])
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout, “main”);

// Define restitution coefficient
T restitutionCoefficient1 = 0.87;
T restitutionCoefficient2 = 0.8;
clout << “Restitution coefficient: ” << restitutionCoefficient1 << std::endl;

// Set bounce to 0
T bounce = 0.0;
clout << “Bounce: ” << bounce << std::endl;

// Calculate cell-wall position
T cellWallPosition = centerY – sphered;
clout << “Cell-wall position: ” << cellWallPosition << std::endl;

// New flag to track bounce
// bool hasBounced = false; // Initially, the particle has not bounced

UnitConverterFromResolutionAndLatticeVelocity<T,DESCRIPTOR> converter(
(int) res, //resolution
( T ) charLatticeVelocity, //charLatticeVelocity
( T ) charPhysLength, //charPhysLength
( T ) charPhysVelocity, //charPhysVelocity
( T ) physViscosity, //physViscosity
( T ) physDensity //physDensity
);
converter.print();

/// === 2rd Step: Prepare Geometry ===
/// Instantiation of a cuboidGeometry with weights
Vector<T,3> origin( 0. );
Vector<T,3> extend( lengthX, lengthY, lengthZ );
IndicatorCuboid3D<T> cuboid(extend, origin);

#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), singleton::mpi().getSize());
#else
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), 7);
#endif
cuboidGeometry.print();

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T,3> superGeometry(cuboidGeometry, loadBalancer, 2);
prepareGeometry(converter, superGeometry);

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

// Prepare lattice
prepareLattice(sLattice, converter, superGeometry);

// Create ParticleSystem
#ifdef PARALLEL_MODE_MPI
SuperParticleSystem<T,PARTICLETYPE> particleSystem(superGeometry);
#else
ParticleSystem<T,PARTICLETYPE> particleSystem;
#endif

//Create particle manager handling coupling, gravity and particle dynamics
ParticleManager<T,DESCRIPTOR,PARTICLETYPE> particleManager(
particleSystem, superGeometry, sLattice, converter, externalAcceleration);

// Create and assign resolved particle dynamics
particleSystem.defineDynamics<
VerletParticleDynamics<T,PARTICLETYPE>>();

// Calculate particle quantities
T epsilon = 0.5*converter.getConversionFactorLength();
Vector<T,3> sphereExtend( sphered );

// Create Particle 1
creators::addResolvedSphere3D(particleSystem, sphereCenter, sphered, epsilon, sphereDensity);
auto particle = particleSystem.get(0);
particle.template setField<MOBILITY,VELOCITY>(sphereVelocity);

// Check ParticleSystem
particleSystem.checkForErrors();

/// === 4th Step: Main Loop with Timer ===
Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel());
timer.start();

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

clout << “MaxIT: ” << converter.getLatticeTime(maxPhysT) << std::endl;
std::ofstream outputFile(“particle_data.dat”);
if (!outputFile.is_open()) {
std::cerr << “Failed to open the file for writing.” << std::endl;
return 1;
}

// Define the total physical time for uniform motion
T uniformMotionTime = 0.02; // 0.02 seconds of uniform motion

// Adjust the time step to ensure that the total motion time is consistent across resolutions
T adjustedTimeStep = uniformMotionTime / (converter.getPhysTime(iTwrite)); // Time adjustment factor based on lattice time step

// Initialize the elapsed time and velocity for uniform motion
T elapsedTime = 0.0;
Vector<T, 3> constantVelocity = sphereVelocity;

size_t writeStep = converter.getLatticeTime(iTwrite);
if (writeStep < 1) writeStep = 1;

T tau = converter.getLatticeRelaxationFrequency();
if (tau <= 0.5) {
clout << “Warning: Tau is too small! Current Tau: ” << tau << std::endl;
exit(1);
}

T cs = 1.0 / sqrt(3.0);
T maxMach = charLatticeVelocity / cs;
if (maxMach >= 0.1) {
clout << “Warning: Ma = ” << maxMach << ” is too large! Reduce charLatticeVelocity.” << std::endl;
exit(1);
}

// Modify the main loop to ensure uniform motion during the specified time
for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT) + 10; ++iT) {

particleManager.execute<
couple_lattice_to_particles<T, DESCRIPTOR, PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
communicate_surface_force<T, PARTICLETYPE>,
#endif
//apply_hydrodynamic_force<T, PARTICLETYPE>,
process_dynamics<T, PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
update_particle_core_distribution<T, PARTICLETYPE>,
#endif
couple_particles_to_lattice<T, DESCRIPTOR, PARTICLETYPE>
>();

//if (iT % converter.getLatticeTime(iTwrite) == 0) {
if (iT % 1 == 0) {
if (singleton::mpi().getRank() == 0) {

elapsedTime = converter.getPhysTime(iT);
if (elapsedTime <= uniformMotionTime) {
auto particle = particleSystem.get(0);
Vector<T, 3> position = particle.template getField<GENERAL, POSITION>();
Vector<T, 3> velocity = constantVelocity;
// Update the particle’s position with the adjusted time step
// position += velocity * adjustedTimeStep; // Adjusting the step for consistent physical time

particle.template setField<GENERAL, POSITION>(position);
particle.template setField<MOBILITY, VELOCITY>(velocity);

elapsedTime += adjustedTimeStep; // Update elapsed time with adjusted step
}

// Execute particle manager
auto particle = particleSystem.get(0); // 从粒子系统中获取粒子
Vector<T, 3> position = particle.template getField<GENERAL, POSITION>();
Vector<T, 3> velocity = particle.template getField<MOBILITY, VELOCITY>();

T bottomPosition = position[1] – sphered;
T epsilon = 2e-4;
if (bottomPosition <= epsilon) {
if (velocity[1] < 0 && bounce < 1) {

std::cout << “Collision detected!” << std::endl;
std::cout << “Velocity before bounce: ” << velocity[1]
<< “, Position before bounce: ” << position[1] << std::endl;

velocity[1] = -velocity[1] * restitutionCoefficient1;
velocity[2] = velocity[2] * restitutionCoefficient2;

// 更新粒子的位置和速度
particle.template setField<GENERAL, POSITION>(position);
particle.template setField<MOBILITY, VELOCITY>(velocity);

bounce += 1;

std::cout << “New velocity after bounce: ” << velocity[1]
<< “, Position after bounce: ” << position[1] << std::endl;
std::cout << “Bounce count: ” << bounce << std::endl;
}
}

if (bounce > 0 && velocity[1] < 0.03) {
velocity[1] *= 2;
particle.template setField<MOBILITY, VELOCITY>(velocity);

}

std::cout << “Time: ” << converter.getPhysTime(iT)
<< ” Current position: ” << position[1] << “, ” << position[2]
<< “, Velocity: ” << velocity[1] << “, ” << velocity[2] << std::endl;

Vector<T, 3> force = particle.template getField<FORCING, FORCE>();
Vector<T, 3> acceleration = particle.template getField<MOBILITY, ACCELERATION_STRD>();

outputFile << converter.getPhysTime(iT) << ” ”
<< position[1] << ” ” << position[2] << ” ”
<< velocity[1] << ” ” << velocity[2] << ” ”
<< acceleration[1] << ” ” << acceleration[2] << ” ”
<< force[1] << ” ” << force[2] << std::endl;
}
}

// Get Results
getResults(sLattice, converter, iT, superGeometry, timer, particleSystem);

// Collide and stream
sLattice.collideAndStreamGpu();

}

outputFile.close();

timer.stop();
timer.printSummary();
}

[2366039968@qq.com]

Hello everyone, I modified the settlingCube3d code to simulate an oblique collision between a particle and a wall. I want to run this modified code on the GPU. My environment has already been configured properly, and the original code can run on the GPU without any issues. However, when I run my modified version in the GPU environment, the GPU memory usage stays at zero. I’m not sure if this is caused by incompatibility in my code.

Below is the code after my modifications.

/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2006-2021 Nicolas Hafen, 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.
*/

/* settlingCube3d.cpp:
* The case examines the settling of a cubical silica particle
* under the influence of gravity.
* The object is surrounded by water in a rectangular domain
* limited by no-slip boundary conditions.
* For the calculation of forces an DNS approach is chosen
* which also leads to a back-coupling of the particle on the fluid,
* inducing a flow.
*
* The simulation is based on the homogenised lattice Boltzmann approach
* (HLBM) introduced in “Particle flow simulations with homogenised
* lattice Boltzmann methods” by Krause et al.
* and extended in “Towards the simulation of arbitrarily shaped 3D particles
* using a homogenised lattice Boltzmann method” by Trunk et al.
* for the simulation of 3D particles.
*
* This example demonstrates the usage of HLBM in the OpenLB framework.
* To improve parallel performance, the particle decomposition scheme
* described in 10.48550/arXiv.2312.14172 is used.
*/

#include “olb3D.h”
#include “olb3D.hh” // use generic version only!

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace olb::particles;
using namespace olb::particles::dynamics;
using namespace std;

using T = FLOATING_POINT_TYPE;

// Define lattice type
typedef PorousParticleD3Q19Descriptor DESCRIPTOR;

// Define particleType
#ifdef PARALLEL_MODE_MPI
// Particle decomposition improves parallel performance
typedef ResolvedDecomposedParticle3D PARTICLETYPE;
#else
typedef ResolvedParticle3D PARTICLETYPE;
#endif

#define WriteVTK
#define WriteGnuPlot

//std::string gnuplotFilename = “gnuplot.dat”;
// Discretization Settings
int res = 800;
T const charLatticeVelocity =0.03;

// Time Settings
T const maxPhysT = 0.15; // max. simulation time in s
T const iTwrite = 0.001; // write out intervall in s

// Domain Settings
T const lengthX = 0.1016;
T const lengthY = 0.1016;
T const lengthZ = 0.1016;

// Fluid Settings
T const physDensity = 1000;
T const physViscosity = 1E-6;
T const sphered = 0.00635/2;
//Particle Settings
T centerX = lengthX*.5;
T centerY =0.0079268;
T centerZ = lengthZ*.5;
T const sphereDensity = 4000;

Vector<T,3> sphereCenter = {centerX,centerY,centerZ};
// External acceleration: gravity
Vector<T,3> externalAcceleration = {.0, .0, -9.81};

// Particle settings: ensure initial velocity is zero
Vector<T,3> sphereVelocity = {0., -0.07884, 0.07884};

// Characteristic Quantities./s
T const charPhysLength = lengthX;
T const charPhysVelocity = 0.2; // Assumed maximal velocity

// Prepare geometry
void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 2);
superGeometry.rename(2, 1, {1, 1, 1});

superGeometry.clean();
superGeometry.innerClean();

superGeometry.checkForErrors();
superGeometry.getStatistics().print();
clout << “Prepare Geometry … OK” << std::endl;
return;
}

// Set up the geometry of the simulation
void prepareLattice(
SuperLatticeGpu<T, DESCRIPTOR> sLattice(superGeometry);
, UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “prepareLattice”); //output information control
clout << “Prepare Lattice …” << std::endl; //out prepare information
clout << “setting Velocity Boundaries …” << std::endl; //velocity boundaries

/// Material=0 –>do nothing
sLattice.defineDynamics<PorousParticleBGKdynamics>(superGeometry, 1); //dynamics module
setBounceBackBoundary(sLattice, superGeometry, 2); // bounce boundary

sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());//relax OMEGA

{
auto& communicator = sLattice.getCommunicator(stage::PostPostProcess());
communicator.requestFields<POROSITY,VELOCITY_NUMERATOR,VELOCITY_DENOMINATOR>();//kongxilv sudufenzi,fenmu
communicator.requestOverlap(sLattice.getOverlap());//overlap area
communicator.exchangeRequests();
}

clout << “Prepare Lattice … OK” << std::endl; //papare ok
}

//Set Boundary Values
void setBoundaryValues(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry)
{
OstreamManager clout(std::cout, “setBoundaryValues”);

if (iT == 0) {
AnalyticalConst3D<T, T> zero(0.);
AnalyticalConst3D<T, T> one(1.);
sLattice.defineField<descriptors::POROSITY>(superGeometry.getMaterialIndicator({0,1,2}), one);
// Set initial condition
AnalyticalConst3D<T, T> ux(0.);
AnalyticalConst3D<T, T> uy(0.);
AnalyticalConst3D<T, T> uz(0.);
AnalyticalConst3D<T, T> rho(1.);
AnalyticalComposed3D<T, T> u(ux, uy, uz);

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU(superGeometry, 1, rho, u);
sLattice.iniEquilibrium(superGeometry, 1, rho, u);

// Make the lattice ready for simulation
sLattice.initialize();
}
}

/// Computes the pressure drop between the voxels before and after the cylinder
void getResults(SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry, Timer<T>& timer,
XParticleSystem<T, PARTICLETYPE>& xParticleSystem )
{
OstreamManager clout(std::cout, “getResults”); //information control

#ifdef WriteVTK
SuperVTMwriter3D<T> vtkWriter(“sedimentation”);
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
SuperLatticePhysExternalPorosity3D<T, DESCRIPTOR> externalPor(sLattice, converter);
vtkWriter.addFunctor(velocity);
vtkWriter.addFunctor(pressure);
vtkWriter.addFunctor(externalPor);

if (iT == 0) {
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}

if (iT % converter.getLatticeTime(iTwrite) == 0) {
vtkWriter.write(iT);
}
#endif

#ifdef PARALLEL_MODE_MPI
communication::forParticlesInSuperParticleSystem<
T, PARTICLETYPE, conditions::valid_particle_centres>(
xParticleSystem,
[&](Particle<T, PARTICLETYPE>& particle,
ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
io::printResolvedParticleInfo(particle);
});
#else
for (std::size_t iP=0; iP<xParticleSystem.size(); ++iP) {
auto particle = xParticleSystem.get(iP);
// io::printResolvedParticleInfo(particle);
}
#endif

}

int main(int argc, char* argv[])
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout, “main”);

// Define restitution coefficient
T restitutionCoefficient1 = 0.87;
T restitutionCoefficient2 = 0.8;
clout << “Restitution coefficient: ” << restitutionCoefficient1 << std::endl;

// Set bounce to 0
T bounce = 0.0;
clout << “Bounce: ” << bounce << std::endl;

// Calculate cell-wall position
T cellWallPosition = centerY – sphered;
clout << “Cell-wall position: ” << cellWallPosition << std::endl;

// New flag to track bounce
// bool hasBounced = false; // Initially, the particle has not bounced

UnitConverterFromResolutionAndLatticeVelocity<T,DESCRIPTOR> converter(
(int) res, //resolution
( T ) charLatticeVelocity, //charLatticeVelocity
( T ) charPhysLength, //charPhysLength
( T ) charPhysVelocity, //charPhysVelocity
( T ) physViscosity, //physViscosity
( T ) physDensity //physDensity
);
converter.print();

/// === 2rd Step: Prepare Geometry ===
/// Instantiation of a cuboidGeometry with weights
Vector<T,3> origin( 0. );
Vector<T,3> extend( lengthX, lengthY, lengthZ );
IndicatorCuboid3D<T> cuboid(extend, origin);

#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), singleton::mpi().getSize());
#else
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getConversionFactorLength(), 7);
#endif
cuboidGeometry.print();

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T,3> superGeometry(cuboidGeometry, loadBalancer, 2);
prepareGeometry(converter, superGeometry);

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

// Prepare lattice
prepareLattice(sLattice, converter, superGeometry);

// Create ParticleSystem
#ifdef PARALLEL_MODE_MPI
SuperParticleSystem<T,PARTICLETYPE> particleSystem(superGeometry);
#else
ParticleSystem<T,PARTICLETYPE> particleSystem;
#endif

//Create particle manager handling coupling, gravity and particle dynamics
ParticleManager<T,DESCRIPTOR,PARTICLETYPE> particleManager(
particleSystem, superGeometry, sLattice, converter, externalAcceleration);

// Create and assign resolved particle dynamics
particleSystem.defineDynamics<
VerletParticleDynamics<T,PARTICLETYPE>>();

// Calculate particle quantities
T epsilon = 0.5*converter.getConversionFactorLength();
Vector<T,3> sphereExtend( sphered );

// Create Particle 1
creators::addResolvedSphere3D(particleSystem, sphereCenter, sphered, epsilon, sphereDensity);
auto particle = particleSystem.get(0);
particle.template setField<MOBILITY,VELOCITY>(sphereVelocity);

// Check ParticleSystem
particleSystem.checkForErrors();

/// === 4th Step: Main Loop with Timer ===
Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel());
timer.start();

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

clout << “MaxIT: ” << converter.getLatticeTime(maxPhysT) << std::endl;
std::ofstream outputFile(“particle_data.dat”);
if (!outputFile.is_open()) {
std::cerr << “Failed to open the file for writing.” << std::endl;
return 1;
}

// Define the total physical time for uniform motion
T uniformMotionTime = 0.02; // 0.02 seconds of uniform motion

// Adjust the time step to ensure that the total motion time is consistent across resolutions
T adjustedTimeStep = uniformMotionTime / (converter.getPhysTime(iTwrite)); // Time adjustment factor based on lattice time step

// Initialize the elapsed time and velocity for uniform motion
T elapsedTime = 0.0;
Vector<T, 3> constantVelocity = sphereVelocity;

size_t writeStep = converter.getLatticeTime(iTwrite);
if (writeStep < 1) writeStep = 1;

T tau = converter.getLatticeRelaxationFrequency();
if (tau <= 0.5) {
clout << “Warning: Tau is too small! Current Tau: ” << tau << std::endl;
exit(1);
}

T cs = 1.0 / sqrt(3.0);
T maxMach = charLatticeVelocity / cs;
if (maxMach >= 0.1) {
clout << “Warning: Ma = ” << maxMach << ” is too large! Reduce charLatticeVelocity.” << std::endl;
exit(1);
}

// Modify the main loop to ensure uniform motion during the specified time
for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT) + 10; ++iT) {

particleManager.execute<
couple_lattice_to_particles<T, DESCRIPTOR, PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
communicate_surface_force<T, PARTICLETYPE>,
#endif
//apply_hydrodynamic_force<T, PARTICLETYPE>,
process_dynamics<T, PARTICLETYPE>,
#ifdef PARALLEL_MODE_MPI
update_particle_core_distribution<T, PARTICLETYPE>,
#endif
couple_particles_to_lattice<T, DESCRIPTOR, PARTICLETYPE>
>();

//if (iT % converter.getLatticeTime(iTwrite) == 0) {
if (iT % 1 == 0) {
if (singleton::mpi().getRank() == 0) {

elapsedTime = converter.getPhysTime(iT);
if (elapsedTime <= uniformMotionTime) {
auto particle = particleSystem.get(0);
Vector<T, 3> position = particle.template getField<GENERAL, POSITION>();
Vector<T, 3> velocity = constantVelocity;
// Update the particle’s position with the adjusted time step
// position += velocity * adjustedTimeStep; // Adjusting the step for consistent physical time

particle.template setField<GENERAL, POSITION>(position);
particle.template setField<MOBILITY, VELOCITY>(velocity);

elapsedTime += adjustedTimeStep; // Update elapsed time with adjusted step
}

// Execute particle manager
auto particle = particleSystem.get(0); // 从粒子系统中获取粒子
Vector<T, 3> position = particle.template getField<GENERAL, POSITION>();
Vector<T, 3> velocity = particle.template getField<MOBILITY, VELOCITY>();

T bottomPosition = position[1] – sphered;
T epsilon = 2e-4;
if (bottomPosition <= epsilon) {
if (velocity[1] < 0 && bounce < 1) {

std::cout << “Collision detected!” << std::endl;
std::cout << “Velocity before bounce: ” << velocity[1]
<< “, Position before bounce: ” << position[1] << std::endl;

velocity[1] = -velocity[1] * restitutionCoefficient1;
velocity[2] = velocity[2] * restitutionCoefficient2;

// 更新粒子的位置和速度
particle.template setField<GENERAL, POSITION>(position);
particle.template setField<MOBILITY, VELOCITY>(velocity);

bounce += 1;

std::cout << “New velocity after bounce: ” << velocity[1]
<< “, Position after bounce: ” << position[1] << std::endl;
std::cout << “Bounce count: ” << bounce << std::endl;
}
}

if (bounce > 0 && velocity[1] < 0.03) {
velocity[1] *= 2;
particle.template setField<MOBILITY, VELOCITY>(velocity);

}

std::cout << “Time: ” << converter.getPhysTime(iT)
<< ” Current position: ” << position[1] << “, ” << position[2]
<< “, Velocity: ” << velocity[1] << “, ” << velocity[2] << std::endl;

Vector<T, 3> force = particle.template getField<FORCING, FORCE>();
Vector<T, 3> acceleration = particle.template getField<MOBILITY, ACCELERATION_STRD>();

outputFile << converter.getPhysTime(iT) << ” ”
<< position[1] << ” ” << position[2] << ” ”
<< velocity[1] << ” ” << velocity[2] << ” ”
<< acceleration[1] << ” ” << acceleration[2] << ” ”
<< force[1] << ” ” << force[2] << std::endl;
}
}

// Get Results
getResults(sLattice, converter, iT, superGeometry, timer, particleSystem);

// Collide and stream
sLattice.collideAndStreamGpu();

}

outputFile.close();

timer.stop();
timer.printSummary();
}

[mathias]

All particle related code parts does not support GPU usage yet.

[Author]

Posts