Problem on Particles-Settling Cube

[Dip]

From the example particles-> settlingcube3d. I converted it to 2d. It works fine, but problem occurs when I tryto modify the code. Like changing the material number to 4 , it doesnt work( I changed all the places equally). I trying it to modify so that the cube is falling through two different fluids, but when I set the lower region to material number 3, after simulation the particle just sat on it. It looks like it only works with material number 1 and 2,2 being the boundary.

#include “olb2D.h”
#include “olb2D.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 T = FLOATING_POINT_TYPE;

// Define lattice type
typedef PorousParticleD2Q9Descriptor DESCRIPTOR;

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

#define WriteVTK

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

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

// Domain Settings
T const lengthX = 0.01;
T const lengthY = 0.01;

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

//Particle Settings
T centerX = lengthX * 0.5;
T centerY = lengthY * 0.5;

T const cubeDensity = 2500;
T const cubeEdgeLength = 0.0025;
Vector<T,2> cubeCenter = {centerX, centerY};
T cubeOrientation = 0.0;
Vector<T,2> cubeVelocity = {0.0, 0.0};
Vector<T,2> externalAcceleration = {0.0, -9.81 * (1.0 – 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,2>& superGeometry) {
OstreamManager clout(std::cout, “prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

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

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,2>& superGeometry) {
OstreamManager clout(std::cout, “prepareLattice”);
clout << “Prepare Lattice …” << std::endl;
clout << “setting Velocity Boundaries …” << std::endl;

// Material=0 –>do nothing
sLattice.defineDynamics<PorousParticleBGKdynamics>(superGeometry, 1);
setBounceBackBoundary(sLattice, 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,2>& superGeometry) {
OstreamManager clout(std::cout, “setBoundaryValues”);

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

// 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,2>& superGeometry, Timer<T>& timer, XParticleSystem<T, PARTICLETYPE>& xParticleSystem) {
OstreamManager clout(std::cout, “getResults”);

#ifdef WriteVTK
SuperVTMwriter2D<T> vtkWriter(“sedimentation”);
SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);
SuperLatticePhysExternalPorosity2D<T, DESCRIPTOR> externalPor(sLattice, converter);
vtkWriter.addFunctor(velocity);
vtkWriter.addFunctor(pressure);
vtkWriter.addFunctor(externalPor);

if (iT == 0) {
// Writes the converter log file
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank2D<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

// 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();

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

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

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T, 2> 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, 2> cubeExtend(cubeEdgeLength);

// Create Particle 1
creators::addResolvedCuboid2D(particleSystem, cubeCenter, cubeExtend, epsilon, cubeDensity, cubeOrientation);

/* // Create Particle 2
cubeCenter = {centerX, lengthY * T(0.51)};
T cubeAngle = 15.; // Set a single angle value for the second particle
creators::addResolvedCuboid2D(particleSystem, cubeCenter, cubeExtend, epsilon, cubeDensity, cubeAngle);*/

// 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();
}

[Adrian]

Did you check the generated material geometry in e.g. Paraview? It is likely that the clean operation will remove the non-material-1 region’s interior by default, you can add further bulk materials as an argument / drop the clean altogether for this simple a geometry.

[Dip]

Yes, after simulation in this code it works fine. But if I change the material number, it is empty.

[Adrian]

Yes, and this case with changed material number is what my comment is a reply to. I repeat: Did you check the geometry and add the new bulk material to the clean method / remove the call altogether?

[Author]

Posts