Commit 28f1e398 authored by Christoph Schwarzmeier's avatar Christoph Schwarzmeier
Browse files

Merge branch 'mr_particle_packing' into 'master'

Showcase for dense particle packing

See merge request walberla/walberla!522
parents 25952342 90801827
......@@ -4,6 +4,7 @@ add_subdirectory( FluidizedBed )
add_subdirectory( HeatConduction )
add_subdirectory( LightRisingParticleInFluidAMR )
add_subdirectory( Mixer )
add_subdirectory( ParticlePacking )
add_subdirectory( PegIntoSphereBed )
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( PhaseFieldAllenCahn )
......
waLBerla_link_files_to_builddir( *.cfg )
if (WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE)
waLBerla_add_executable ( NAME ParticlePacking
FILES ParticlePacking.cpp
DEPENDS blockforest core mesa_pd postprocessing sqlite vtk )
endif()
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla 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 3 of
// the License, or (at your option) any later version.
//
// waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file DiameterDistribution.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include <random>
#include "Utility.h"
namespace walberla {
namespace mesa_pd {
class DiameterGenerator
{
public:
virtual real_t get() = 0;
virtual ~DiameterGenerator() = default;
};
class LogNormal: public DiameterGenerator
{
public:
LogNormal(real_t mu, real_t var, uint_t seed = 42)
: gen_(seed)
{
auto m = std::log(mu*mu/std::sqrt(var+mu*mu));
auto s = std::sqrt(log(var/(mu*mu) + 1_r));
dist_ = std::lognormal_distribution<> (m, s);
}
real_t get() override
{
return real_c(dist_(gen_));
}
private:
std::mt19937 gen_;
std::lognormal_distribution<> dist_;
};
class Uniform: public DiameterGenerator
{
public:
Uniform(real_t diameter)
: diameter_(diameter)
{ }
real_t get() override
{
return diameter_;
}
private:
real_t diameter_;
};
class DiscreteSieving: public DiameterGenerator
{
public:
DiscreteSieving(std::vector<real_t> diameters, std::vector<real_t> massFractions, uint_t seed,
real_t normalParticleVolume, real_t totalParticleMass = 1_r, real_t particleDensity = 1_r)
: diameters_(diameters), gen_(seed)
{
WALBERLA_CHECK_EQUAL(diameters.size(), massFractions.size(), "Number of entries in diameter and mass-fraction array has to be the same!");
WALBERLA_CHECK_FLOAT_EQUAL(real_t(1), std::accumulate(massFractions.begin(), massFractions.end(), real_t(0)), "Sum of mass fractions has to be 1!");
auto particleNumbers = transferMassFractionsToParticleNumbers(massFractions, diameters, normalParticleVolume, totalParticleMass, particleDensity);
std::string outString = "Discrete Sieving: Expected particle numbers per diameter: | ";
bool issueWarning = false;
for(const auto & p : particleNumbers){
issueWarning |= (p > 0_r && p < 1_r);
outString += std::to_string(int(p)) + " | ";
}
auto totalParticles = std::accumulate(particleNumbers.begin(), particleNumbers.end(), real_t(0));
outString += " -> total = " + std::to_string(int(totalParticles));
WALBERLA_LOG_INFO_ON_ROOT(outString);
if(issueWarning) WALBERLA_LOG_WARNING_ON_ROOT("Attention: The simulated particle mass is not enough to have at least a single particle of all size classes!");
dist_ = std::discrete_distribution<uint_t>(particleNumbers.begin(), particleNumbers.end());
}
real_t get() override
{
return diameters_[dist_(gen_)];
}
private:
std::vector<real_t> diameters_;
std::mt19937 gen_;
std::discrete_distribution<uint_t> dist_;
};
class ContinuousSieving: public DiameterGenerator
{
public:
ContinuousSieving(std::vector<real_t> sieveSizes, std::vector<real_t> massFractions, uint_t seed,
real_t normalParticleVolume, real_t totalParticleMass = 1_r, real_t particleDensity = 1_r)
: gen_(seed)
{
WALBERLA_CHECK_EQUAL(sieveSizes.size(), massFractions.size()+1, "Number of entries in sieves has to be one larger than the mass-fraction array!");
WALBERLA_CHECK_FLOAT_EQUAL(real_t(1), std::accumulate(massFractions.begin(), massFractions.end(), real_t(0)), "Sum of mass fractions has to be 1!");
auto meanDiameters = getMeanDiametersFromSieveSizes(sieveSizes);
auto particleNumbers = transferMassFractionsToParticleNumbers(massFractions, meanDiameters, normalParticleVolume, totalParticleMass, particleDensity);
std::string outString = "Continuous Sieving: Expected particle numbers per diameter: | ";
bool issueWarning = false;
for(const auto & p : particleNumbers){
issueWarning |= (p > 0_r && p < 1_r);
outString += std::to_string(int(p)) + " | ";
}
auto totalParticles = std::accumulate(particleNumbers.begin(), particleNumbers.end(), real_t(0));
outString += " -> total = " + std::to_string(int(totalParticles));
WALBERLA_LOG_INFO_ON_ROOT(outString);
if(issueWarning) WALBERLA_LOG_WARNING_ON_ROOT("Attention: The simulated particle mass is not enough to have at least a single particle of all size classes!");
std::vector<real_t> logSieveSizes(sieveSizes.size(),0_r); // = phi-scale -> uniformly distributed
for(uint_t i = 0; i < logSieveSizes.size(); ++i)
{
logSieveSizes[i] = - std::log2(sieveSizes[i]);
}
dist_ = std::piecewise_constant_distribution<real_t>(logSieveSizes.begin(), logSieveSizes.end(), particleNumbers.begin());
}
real_t get() override
{
return std::pow(2_r,-dist_(gen_)); //re-transfer from phi-scale to actual diameter
}
private:
std::mt19937 gen_;
std::piecewise_constant_distribution<real_t> dist_;
};
} // namespace mesa_pd
} // namespace walberla
This diff is collapsed.
ParticlePacking
{
domainSetup periodic; // container, periodic
domainWidth 0.104; // m
domainHeight 1; // m
particleDensity 2650; // kg/m^3, 2500 (spheres), 2650 (real sediments);
ambientDensity 1000; // kg/m^3
gravitationalAcceleration 9.81; // m/s^2
particleDistribution SievingCurve; // size distribution, see 'Distribution' block for options
particleShape Mesh; // see 'Shape' block
limitVelocity -1; // m/s, negative switches limiting off
initialVelocity 1; // m/s
initialGenerationHeightRatioStart 0.1; // -
initialGenerationHeightRatioEnd 1; // -
generationSpacing 0.010; // m
scaleGenerationSpacingWithForm true;
generationHeightRatioStart 0.6; // -
generationHeightRatioEnd 1; // -
totalParticleMass 4; // kg
numBlocksPerDirection <3,3,1>;
visSpacing 0.01; // s
infoSpacing 0.01; // s
loggingSpacing 0.01; // s
terminalVelocity 1e-3; // m/s
terminalRelativeHeightChange 1e-5; // -
terminationCheckingSpacing 0.01; // s
minimalTerminalRunTime .0; // s, minimal runtime after generation & shaking
solver DEM; // particle simulation approach, see 'Solver' block for options
shaking true; // see 'Shaking' block
velocityDampingCoefficient 1e-3; // continuous reduction of velocity in last simulation phase
useHashGrids false;
particleSortingSpacing 1000; // time steps, non-positive values switch sorting off, performance optimization
}
Shaking
{
amplitude 3e-4; // m
period 0.025; // s
duration 6.0; // s, duration of shaking AFTER creation of all particles
activeFromBeginning true;
}
Solver
{
coefficientOfRestitution 0.1; // -
frictionCoefficient 0.5; // -
dt 5e-6; // s, time step size
DEM
{
collisionTime 20e-5; // s
poissonsRatio 0.22; // -
}
HCSITS
{
errorReductionParameter 0.8;
relaxationParameter 0.75;
numberOfIterations 10;
relaxationModel InelasticGeneralizedMaximumDissipationContact;
}
}
Distribution
{
randomSeed 41; // if negative, seed is randomly chosen
Uniform{
diameter 11e-3; // m
}
LogNormal{
mu 1e-3; // m
variance 1e-7;
}
DiameterMassFractions
{
// from Schruff, 2016
//diameters 1.5e-3 2e-3 3e-3 4e-3 6e-3 8e-3 11e-3 16e-3 22e-3;
//massFractions .00 .00 1.00 .00 .00 .00 .00 .00 .00; // I, U
//massFractions .00 .00 .00 1.00 .00 .00 .00 .00 .00; // II, U
//massFractions .00 .00 .00 .00 1.00 .00 .00 .00 .00; // III, U
//massFractions .00 .00 .50 .00 .50 .00 .00 .00 .00; // IV, B
//massFractions .00 .00 .00 .50 .00 .50 .00 .00 .00; // V, B
//massFractions .00 .00 .00 .00 .00 .50 .00 .50 .00; // VI, B
//massFractions .00 .00 .25 .50 .25 .00 .00 .00 .00; // VII, T
//massFractions .00 .00 .00 .25 .50 .25 .00 .00 .00; // VIII, T
//massFractions .00 .00 .00 .00 .25 .50 .25 .00 .00; // IX, T
//massFractions .025 .050 .100 .200 .250 .200 .100 .050 .025; // X, LN
// Liang, 2015, discrete spheres (3,4,6,8,11,16,22 mm)
diameters 3e-3 4e-3 6e-3 8e-3 11e-3 16e-3 22e-3;
//massFractions 0 0 1 0 0 0 0; // A
//massFractions .0 .0 .21 .58 .21 0 0; // B
//massFractions .0 .06 .24 .4 .24 .06 .0; // C
//massFractions .04 .11 .22 .26 .22 .11 .04; // D
//massFractions .08 .13 .08 .06 .18 .29 .18; // E
massFractions .13 .21 .13 .06 .13 .21 .13; // F
//massFractions .18 .29 .18 .06 .08 .13 .08; // G
}
SievingCurve
{
// Frings, 2011, Verification, slightly truncated, table 2
//sieveSizes 0.063e-3 0.09e-3 0.125e-3 0.18e-3 0.25e-3 0.355e-3 0.5e-3 0.71e-3 1e-3 1.4e-3 2e-3 2.8e-3 4e-3 5.6e-3 8e-3 11.2e-3 16e-3 22.4e-3 31.5e-3 45e-3 63e-3;
//massFractions .00 .00 .01 .01 .05 .15 .12 .06 .03 .03 .03 .05 .06 .07 .07 .07 .06 .06 .04 .03; // case 45
// Lian, 2015, DigiPack, values from porosity report (Versuch3), in cylinder with diam 104mm
sieveSizes 2.8e-3 4e-3 5.6e-3 8e-3 11.2e-3 16e-3 22.4e-3 31.5e-3;
//massFractions 0 0 1 0 0 0 0; // A
//massFractions .0 .0 .21 .58 .21 0 0; // B
//massFractions .0 .06 .24 .4 .24 .06 .0; // C
//massFractions .04 .11 .22 .26 .22 .11 .04; // D
//massFractions .08 .13 .08 .06 .18 .29 .18; // E
massFractions .13 .21 .13 .06 .13 .21 .13; // F
//massFractions .18 .29 .18 .06 .08 .13 .08; // G
//massFractions 0 0 0 0 0 0 1; // D
useDiscreteForm false; // only use average sieving diameters
}
}
Shape
{
scaleMode sieveLike; // sphereEquivalent, sieveLike
Sphere
{
}
Ellipsoid
{
semiAxes <50,5,1>; // will be scaled to obtain desired size
}
EquivalentEllipsoid
{
// specify either a mesh file or a folder containing mesh files ( have to be .obj files)
//path example_meshes/sediment_scan_0.obj;
path example_meshes;
}
EllipsoidFormDistribution
{
elongationMean 0.4;
elongationStdDev 0.2;
flatnessMean 0.4;
flatnessStdDev 0.2;
}
Mesh
{
// specify either a mesh file or a folder containing mesh files ( have to be .obj files)
//path example_meshes/sediment_scan_0.obj;
path example_meshes;
}
MeshFormDistribution
{
path example_meshes; // meshes used as basis, are then scaled to match size and form
elongationMean 0.4; //0.682;
elongationStdDev 0.2; //0.139;
flatnessMean 0.4; //0.65;
flatnessStdDev 0.2; //0.175;
}
UnscaledMeshesPerFraction
{
// expects subfolders 0, 1,... in given folder that contain a set of meshes
// those meshes are taken as is and not scaled during creation
// the mass fraction from Distribution/SievingCurve is used to determine the generation probabilities per fraction
folder meshes_collection_unscaled;
}
}
Evaluation
{
vtkFolder vtk_out_container;
vtkFinalFolder vtk_files;
//Rhein, Frings, 2011, Verification, Table 2
histogramBins 200e-3 125e-3 90e-3 63e-3 45e-3 31.5e-3 22.4e-3 16e-3 11.2e-3 8e-3 5.6e-3 4e-3 2.8e-3 2e-3 1.4e-3 1e-3 0.71e-3 0.5e-3 0.355e-3 0.25e-3 0.18e-3 0.125e-3 0.09e-3 0.063e-3;
//for horizontal layer evaluation (only spheres)
porosityProfileFolder porosity_profiles;
layerHeight 1e-3; // m
sqlDBFileName db_ParticlePacking.sqlite;
}
This diff is collapsed.
This diff is collapsed.
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla 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 3 of
// the License, or (at your option) any later version.
//
// waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file Utility.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "mesa_pd/data/ParticleStorage.h"
#include <iterator>
#include <algorithm>
#include <functional>
namespace walberla {
namespace mesa_pd {
class ParticleAccessorWithShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& particleStorage)
: ParticleAccessor(particleStorage)
{}
walberla::real_t const & getInvMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvMass();}
Mat3 const & getInvInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvInertiaBF();}
Mat3 const & getInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInertiaBF();}
data::BaseShape* getShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx).get();}
};
template< typename T>
std::vector<T> parseStringToVector(std::string inputStr)
{
std::istringstream iss(inputStr);
return std::vector<T>{std::istream_iterator<T>(iss),std::istream_iterator<T>()};
}
real_t radiusFromSphereVolume(real_t volume)
{
return std::cbrt(real_t(3) / ( real_t(4) * math::pi) * volume);
}
real_t diameterFromSphereVolume(real_t volume)
{
return real_t(2) * radiusFromSphereVolume(volume);
}
uint_t getIndexOfSecondSemiAxis(Vector3<real_t> semiAxes)
{
std::set<uint_t> indices = {0,1,2};
indices.erase(semiAxes.indexOfMin());
indices.erase(semiAxes.indexOfMax());
return *indices.begin();
}
void sortVector(Vector3<real_t> & v)
{
std::sort(v.data(),v.data()+3);
}
real_t sizeFromSemiAxes(Vec3 semiAxes)
{
sortVector(semiAxes);
return 2_r * std::sqrt((semiAxes[0] * semiAxes[0] + semiAxes[1] * semiAxes[1]) / 2_r); // factor of 2 because those are SEMI axis
}
// assuming a mass-equivalent ellipsoid
Vec3 semiAxesFromInertiaTensor(const Matrix3<real_t> & inertiaTensor, real_t mass)
{
Vec3 semiAxes;
semiAxes[0] = std::sqrt( (- inertiaTensor(0,0) + inertiaTensor(1,1) + inertiaTensor(2,2)) / (2_r * mass / 5_r) );
semiAxes[1] = std::sqrt( (- inertiaTensor(1,1) + inertiaTensor(2,2) + inertiaTensor(0,0)) / (2_r * mass / 5_r) );
semiAxes[2] = std::sqrt( (- inertiaTensor(2,2) + inertiaTensor(0,0) + inertiaTensor(1,1)) / (2_r * mass / 5_r) );
return semiAxes;
}
std::vector<real_t> getMeanDiametersFromSieveSizes(std::vector<real_t> sieveSizes)
{
//since grain sizes are logarithmically distributed, it is practice to use the geometric mean
std::vector<real_t> meanDiameters(sieveSizes.size()-1,0_r);
for(uint_t i = 0; i < meanDiameters.size(); ++i)
{
meanDiameters[i] = std::sqrt(sieveSizes[i] * sieveSizes[i+1]);
}
return meanDiameters;
}
// if totalMass and density are given actual numbers, the resulting particle numbers are a good estimate for the actual numbers
// else, the resulting numbers are directly proportional to the actual ones, which is sufficient to define the distributions
std::vector<real_t> transferMassFractionsToParticleNumbersFromAvgVolumes(std::vector<real_t> massFractions, std::vector<real_t> avgVolumePerSizeFraction, real_t totalMass = real_t(1), real_t density = real_t(1) )
{
WALBERLA_CHECK_EQUAL(avgVolumePerSizeFraction.size(), massFractions.size(), "Number of entries in volume and mass-fraction array has to be the same!");
std::vector<real_t> particleNumbers(massFractions.size(), real_t(0));
for(uint_t n = 0; n < massFractions.size(); ++n )
{
if(avgVolumePerSizeFraction[n] > 0_r) particleNumbers[n] = totalMass * massFractions[n] / (density * avgVolumePerSizeFraction[n]);
}
return particleNumbers;
}
// note: normalVolume is the volume of a typical particle with diameter = 1. For sphere: PI / 6
std::vector<real_t> transferMassFractionsToParticleNumbers(std::vector<real_t> massFractions, std::vector<real_t> diameters, real_t normalVolume = math::pi / real_t(6), real_t totalMass = real_t(1), real_t density = real_t(1) )
{
WALBERLA_CHECK_EQUAL(diameters.size(), massFractions.size(), "Number of entries in diameter and mass-fraction array has to be the same!");
std::vector<real_t> avgVolumePerSizeFraction(massFractions.size(), real_t(0));
for(uint_t n = 0; n < massFractions.size(); ++n )
{
avgVolumePerSizeFraction[n] = normalVolume * diameters[n] * diameters[n] * diameters[n];
}
return transferMassFractionsToParticleNumbersFromAvgVolumes(massFractions, avgVolumePerSizeFraction, totalMass, density);
}
real_t computePercentileFromSieveDistribution(std::vector<real_t> diameters, std::vector<real_t> massFractions, real_t percentile)
{
WALBERLA_CHECK(percentile > 0_r && percentile < 100_r, "Percentile is a value between 0 and 100");
if(diameters[0] > diameters[1])
{
// reverse order to have it ascending
std::reverse(diameters.begin(), diameters.end());
std::reverse(massFractions.begin(), massFractions.end());
}
std::vector<real_t> cdf(massFractions.size(),0_r);
std::partial_sum(massFractions.begin(), massFractions.end(), cdf.begin());
for(uint_t i = 0; i < cdf.size()-1; ++i)
{
if(cdf[i] <= percentile/100_r && percentile/100_r <= cdf[i+1] )
{
real_t f_small = cdf[i];
real_t f_large = cdf[i+1];
if(f_small <= 0.000001_r && f_large>= 0.999999_r)
{
// special case of uniform distribution -> percentile is this value
return diameters[i+1];
}
real_t phi_small = - std::log2(diameters[i]);
real_t phi_large = - std::log2(diameters[i+1]);
// logarithmic interpolation of diameter value
real_t phi_percentile = phi_small + (phi_large - phi_small) / (f_large - f_small) * (percentile / 100_r - f_small);
return std::pow(2_r,-phi_percentile);
}
}
return diameters[0];
}
auto createPlane( std::shared_ptr<data::ParticleStorage> particleStorage,
const Vec3& pos,
const Vec3& normal)
{
auto p = particleStorage->create(true);
p->setPosition( pos );
p->setBaseShape( std::make_shared<data::HalfSpace>( normal ) );
p->getBaseShapeRef()->updateMassAndInertia(real_t(1));
p->setOwner( walberla::mpi::MPIManager::instance()->rank() );
p->setType( 0 );
p->setInteractionRadius(std::numeric_limits<real_t>::infinity());
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::GLOBAL);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::INFINITE);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::FIXED);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::NON_COMMUNICATING);
return p;
}
auto createCylindricalBoundary( std::shared_ptr<data::ParticleStorage> particleStorage,
const Vec3& pos,
const Vec3& axis, real_t radius)
{
auto p = particleStorage->create(true);
p->setPosition( pos );
p->setBaseShape( std::make_shared<data::CylindricalBoundary>( radius, axis ) );
p->getBaseShapeRef()->updateMassAndInertia(real_t(1));
p->setOwner( walberla::mpi::MPIManager::instance()->rank() );
p->setType( 0 );
p->setInteractionRadius(std::numeric_limits<real_t>::infinity());
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::GLOBAL);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::INFINITE);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::FIXED);
data::particle_flags::set(p->getFlagsRef(), data::particle_flags::NON_COMMUNICATING);
return p;
}
void writeParticleInformationToFile(const std::string& filename, const std::string& particleInfoStr, bool logToProcessLocalFiles)
{
if(logToProcessLocalFiles)
{
std::ofstream file;
file.open( filename.c_str());
file << particleInfoStr;
file.close();
}else{
walberla::mpi::writeMPITextFile( filename, particleInfoStr );
}