PSD Class Reference

Contains a vector with radii and probabilities of a user defined particle size distribution (PSD) More...

#include <PSD.h>

Public Types
enum class	TYPE {   CUMULATIVE_NUMBER_DISTRIBUTION , CUMULATIVE_LENGTH_DISTRIBUTION , CUMULATIVE_VOLUME_DISTRIBUTION , CUMULATIVE_AREA_DISTRIBUTION ,   PROBABILITYDENSITY_NUMBER_DISTRIBUTION , PROBABILITYDENSITY_LENGTH_DISTRIBUTION , PROBABILITYDENSITY_AREA_DISTRIBUTION , PROBABILITYDENSITY_VOLUME_DISTRIBUTION }
	Enum class which stores the possible types of CDFs and PDFs. Particle size distributions can be represented by different probabilities based on number of particles, length of particles, surface area of particles or volume of particles. More...

Public Member Functions
	PSD ()
	Constructor; sets everything to 0 or default. More...
	PSD (std::unique_ptr< PSD > psdPtr)
	Constructor; sets everything to 0 or default. This constructor is mainly used in MercuryPBM where smart pointers are used. More...
	PSD (const PSD &other)
	Copy constructor with deep copy. More...
	~PSD ()
	Destructor; default destructor. More...
PSD *	copy () const
	Creates a copy on the heap and returns a pointer. More...
void	printPSD () const
	Prints radii and probabilities of the PSD vector. More...
Mdouble	drawSample ()
	Draw a sample radius from a CUMULATIVE_NUMBER_DISTRIBUTION. More...
Mdouble	insertManuallyByVolume (Mdouble volume)
	Draw sample radius manually per size class and check the volumeAllowed of each size class to insert the PSD as accurate as possible. More...
void	validateCumulativeDistribution ()
	Validates if a CDF starts with zero and adds up to unity. More...
void	validateProbabilityDensityDistribution ()
	Validates if the integral of the PDF equals to unity. More...
void	setPSDFromVector (std::vector< DistributionElements > psd, TYPE PSDType)
	Sets the PSD from a vector with DistributionElements. More...
void	setPSDFromCSV (const std::string &fileName, TYPE PSDType, bool headings=false, Mdouble unitScalingFactorRadii=1.0)
	read in the PSD vector with probabilities and radii saved in a .csv file. More...
void	setDistributionUniform (Mdouble radMin, Mdouble radMax, int numberOfBins)
	create a PSD vector for a uniform distribution. More...
void	setDistributionNormal (Mdouble mean, Mdouble standardDeviation, int numberOfBins)
	create a PSD vector for a normal distribution. More...
void	setDistributionLogNormal (Mdouble mean, Mdouble standardDeviation, int numberOfBins)
	create a PSD vector for a normal distribution. More...
void	setDistributionPhiNormal (Mdouble D50, Mdouble standardDeviationinPhi, int numberOfBins)
	create a PSD vector for a normal distribution in Phi Units that has the demanded D50. More...
void	convertProbabilityDensityToCumulative ()
	Converts a PDF to a CDF by integration. More...
void	convertCumulativeToProbabilityDensity ()
	Converts a CDF to a PDF by derivation. More...
void	convertProbabilityDensityToProbabilityDensityNumberDistribution (TYPE PDFType)
	Convert any PDF to a PROBABILITYDENSITY_NUMBER_DISTRIBUTION. More...
void	convertProbabilityDensityNumberDistributionToProbabilityDensity (TYPE PDFType)
	Convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to any PDF. More...
MERCURYDPM_DEPRECATED void	convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution ()
	convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION. More...
void	convertCumulativeToCumulativeNumberDistribution (TYPE CDFType)
	Convert any CDF to a CUMULATIVE_NUMBER_DISTRIBUTION. More...
void	convertCumulativeNumberDistributionToCumulative (TYPE CDFType)
	Convert a CUMULATIVE_NUMBER_DISTRIBUTION to any CDF. More...
MERCURYDPM_DEPRECATED void	cutoffCumulativeNumber (Mdouble quantileMin, Mdouble quantileMax, Mdouble minPolydispersity=0.1)
	cutoff the PSD at given quantiles. More...
void	cutoffCumulativeDistribution (TYPE CDFType, Mdouble quantileMin, Mdouble quantileMax, Mdouble minPolydispersity=0.1)
	Cut off the PSD at the given quantiles, with the PSD in the provided CDF form. More...
MERCURYDPM_DEPRECATED void	cutoffAndSqueezeCumulative (Mdouble quantileMin, Mdouble quantileMax, Mdouble squeeze, Mdouble minPolydispersity=0.1)
	cutoff the PSD at given quantiles and make it less polydisperse by squeezing it. More...
void	cutoffAndSqueezeCumulativeDistribution (TYPE CDFType, Mdouble quantileMin, Mdouble quantileMax, Mdouble squeeze, Mdouble minPolydispersity=0.1)
	Cut off the PSD at the given quantiles, with the PSD in the provided CDF form, and make it less polydisperese by squeezing it. More...
Mdouble	getMinRadius () const
	Get smallest radius of the PSD. More...
Mdouble	getMaxRadius () const
	Get largest radius of the PSD. More...
std::vector< DistributionElements >	getParticleSizeDistribution () const
	Get the PSD vector. More...
std::vector< DistributionElements >	getParticleSizeDistributionByType (TYPE psdType, Mdouble scalingFactor=1.0) const
	Gets the PSD vector, converted to the preferred type. More...
MERCURYDPM_DEPRECATED void	setParticleSizeDistribution (std::vector< DistributionElements >)
	set the PSD by a suitable vector. More...
void	scaleParticleSize (double scale)
	Scales all particle sizes by a factor. More...
double	scaleParticleSizeAuto (int numberOfParticles, double targetVolume, bool allowScaleDown=false)
	Scales all particle sizes, such that the total volume of N particles approximately equals the target volume. More...
int	getInsertedParticleNumber () const
	Get the number of particles already inserted into the simulation. More...
void	decrementNParticlesPerClass ()
	Decrement nParticlesPerClass_ counter. More...
void	decrementVolumePerClass (Mdouble volume)
	Decrement volumePerClass_ counter. More...
Mdouble	getRadius (int index) const
	get a radius at a certain index of the particleSizeDistribution_ More...
Mdouble	getNumberDx (Mdouble x) const
	Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the number based PSD. More...
Mdouble	getLengthDx (Mdouble x) const
	Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the length based PSD. More...
Mdouble	getAreaDx (Mdouble x) const
	Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the area based PSD. More...
Mdouble	getVolumeDx (Mdouble x) const
	Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the volume based PSD. More...
Mdouble	getRadiusByQuantile (Mdouble quantile) const
	Calculate the quantile of the PSD. More...
Mdouble	getQuantileByRadius (Mdouble radius) const
	Calculates the quantile corresponding to a certain radius. More...
Mdouble	getVolumetricMeanRadius () const
	get a volumetric mean radius of the PSD. More...
Mdouble	getSizeRatio () const
	get the size ratio (width) of the PSD. More...
void	cutHighSizeRatio ()
	Check if the size ratio is too high and cut it. More...
std::vector< std::pair< Mdouble, Mdouble > >	convertPdfPhiToPdfMeter (Mdouble meaninPhi, Mdouble standardDeviationinPhi, int numberOfBins)
	compute central momenta of the user defined PSD. ‍/ void computeCentralMomenta(); More...
Mdouble	computeD50 (std::vector< std::pair< Mdouble, Mdouble > > vectorOfPair)
	compute the mass-based D_50 of the paired diameter and probabilityDensity. More...
void	setFixedSeed (int seed)
	set a fixed seed for the random number generator; this is used for the PSDSelfTest to reproduce results More...

Static Public Member Functions
static PSD	getDistributionNormal (Mdouble mean, Mdouble standardDeviation, int numberOfBins)
static PSD	getDistributionLogNormal (Mdouble mean, Mdouble standardDeviation, int numberOfBins)
static std::vector< Mdouble >	linspace (Mdouble Min, Mdouble Max, int numberOfBins)
	create a vector of linearly spaced values. More...

Private Attributes
std::vector< DistributionElements >	particleSizeDistribution_
std::vector< int >	nParticlesPerClass_
std::vector< Mdouble >	volumePerClass_
int	index_
RNG	random_

Friends
bool	operator< (const DistributionElements &l, const DistributionElements &r)
	determines if a certain value of the PSD vector is lower than another one. Used for std::lower_bound() More...
bool	operator< (const DistributionElements &l, Mdouble r)
	determines if a certain value of the PSD vector is lower than a double. More...
std::ostream &	operator<< (std::ostream &os, DistributionElements &p)
	Writes to output stream. More...
std::istream &	operator>> (std::istream &is, DistributionElements &p)
	Reads from input stream. More...
Mdouble	operator== (DistributionElements l, Mdouble r)
	Determines if a certain value of the PSD vector is equal to a double. More...

Detailed Description

Contains a vector with radii and probabilities of a user defined particle size distribution (PSD)

Stores radii and probabilities of a particle size distribution (PSD) in a vector of type DistributionElements and converts them to a PSD which can be used by insertionBoundaries which insert particles into a simulation:

Cumulative distribution function (CDF): gives percentage of particles p_i whose radius r is less than a certain radius r_i. It requires p_0 = 0 and p_end = 1. The CDF's probabilities can be number, length, area or volume based. The cumulative number distribution function (CUMULATIVE_NUMBER_DISTRIBUTION) is used as PSD in MercuryDPM, as it can be interpreted as the probability that a particle's radius is less than r_i: CDF(r<r_i) = p_i.

Probability density function (PDF): p_i is the percentage of particles whose radius is between r_i-1 and r_i. It requires p_0=0 and sum(p_i)=1. The PDF's probabilities can also be number, length, area or volume based. PDF's are utilized to convert any type of PDF to a probability number density function (PROBABILITYDENSITY_NUMBER_DISTRIBUTION) which in a next step are converted to the default CUMULATIVE_NUMBER_DISTRIBUTION.

Sieve data: p_i is the percentage of particles whose radius is between r_i and r_i+1. It requires sum(p_i)=1. relation to PDF: p_i = p(PDF)_i-1. Sieve data is not yet used in this class.

Default distribution is the CUMULATIVE_NUMBER_DISTRIBUTION.

Member Enumeration Documentation

TYPE

enum PSD::TYPE

strong

Enum class which stores the possible types of CDFs and PDFs. Particle size distributions can be represented by different probabilities based on number of particles, length of particles, surface area of particles or volume of particles.

Enumerator
CUMULATIVE_NUMBER_DISTRIBUTION
CUMULATIVE_LENGTH_DISTRIBUTION
CUMULATIVE_VOLUME_DISTRIBUTION
CUMULATIVE_AREA_DISTRIBUTION
PROBABILITYDENSITY_NUMBER_DISTRIBUTION
PROBABILITYDENSITY_LENGTH_DISTRIBUTION
PROBABILITYDENSITY_AREA_DISTRIBUTION
PROBABILITYDENSITY_VOLUME_DISTRIBUTION

    {
        CUMULATIVE_NUMBER_DISTRIBUTION,
        CUMULATIVE_LENGTH_DISTRIBUTION,
        CUMULATIVE_VOLUME_DISTRIBUTION,
        CUMULATIVE_AREA_DISTRIBUTION,
        PROBABILITYDENSITY_NUMBER_DISTRIBUTION,
        PROBABILITYDENSITY_LENGTH_DISTRIBUTION,
        PROBABILITYDENSITY_AREA_DISTRIBUTION,
        PROBABILITYDENSITY_VOLUME_DISTRIBUTION
    };

Constructor & Destructor Documentation

PSD() [1/3]

PSD::PSD

(

)

Constructor; sets everything to 0 or default.

Default constructor; sets every data member to 0 or default.

{
//    for (auto& momenta: momenta_)
//        momenta = 0;
    index_ = 0;
    // default seed is zero for reproducibility
    random_.setRandomSeed(0);
}

References index_, random_, and RNG::setRandomSeed().

Referenced by copy().

PSD() [2/3]

PSD::PSD ( std::unique_ptr< PSD >  psdPtr )

inline

Constructor; sets everything to 0 or default. This constructor is mainly used in MercuryPBM where smart pointers are used.

Parameters

[in]

psdPtr

smart pointer to a PSD object.

                                   : particleSizeDistribution_(psdPtr->particleSizeDistribution_),
                                       nParticlesPerClass_(psdPtr->nParticlesPerClass_),
                                       volumePerClass_(psdPtr->volumePerClass_),
                                       index_(psdPtr->index_),
                                       random_(psdPtr->random_) {}

PSD() [3/3]

PSD::PSD

(

const PSD &

other

)

Copy constructor with deep copy.

Copy constructor

{
    particleSizeDistribution_ = other.particleSizeDistribution_;
//    momenta_ = other.momenta_;
    index_ = other.index_;
    random_ = other.random_;
}

References index_, particleSizeDistribution_, and random_.

~PSD()

PSD::~PSD ( )

default

Destructor; default destructor.

Destructor. Since there are no pointers in this class, there is no need for any actions here.

Member Function Documentation

computeD50()

Mdouble PSD::computeD50

(

std::vector< std::pair< Mdouble, Mdouble > >

vectorOfPair

)

compute the mass-based D_50 of the paired diameter and probabilityDensity.

computes the mass-based D50 of the paired diameter and probabilityDensity

Parameters

[in]

vectorOfPair

vector of pair containing paired data of diameter and probabilityDensity

    {
        Mdouble D_50 = 0.0;
        int numberOfBins = vectorOfPair.size();
        std::vector<Mdouble> Vcum(numberOfBins, 0);
        Mdouble Vint=0;
        for (int i = 0; i < Vcum.size()-1; i++)
        {
                Mdouble dD = fabs(vectorOfPair[i+1].first-vectorOfPair[i].first);
                Mdouble V = 0.5*(pow(vectorOfPair[i].first, 3.0) * vectorOfPair[i].second +
                        pow(vectorOfPair[i+1].first, 3.0) * vectorOfPair[i+1].second) * dD;
                Vint = Vint + V;
                Vcum[i+1] = Vint;
        }
 
        for(int j=0; j<Vcum.size()-1; ++j){
            if(Vcum[j]/Vint*100<=50 && Vcum[j+1]/Vint*100>50){
                D_50 = vectorOfPair[j].first + fabs(vectorOfPair[j+1].first-vectorOfPair[j].first)/(Vcum[j+1]/Vint*100-Vcum[j]/Vint*100)
                        *(50-Vcum[j]/Vint*100);
            }
        }
        return D_50;
    }

References boost::multiprecision::fabs(), i, j, Eigen::bfloat16_impl::pow(), and V.

Referenced by setDistributionPhiNormal().

convertCumulativeNumberDistributionToCumulative()

void PSD::convertCumulativeNumberDistributionToCumulative

(

TYPE

CDFType

)

Convert a CUMULATIVE_NUMBER_DISTRIBUTION to any CDF.

Parameters

CDFType

The CDF type to convert to.

{
    // Ignore the default case.
    if (CDFType == TYPE::CUMULATIVE_NUMBER_DISTRIBUTION)
        return;
 
    convertCumulativeToProbabilityDensity();
 
    if (CDFType == TYPE::CUMULATIVE_LENGTH_DISTRIBUTION)
        convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION);
    else if (CDFType == TYPE::CUMULATIVE_AREA_DISTRIBUTION)
        convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION);
    else if (CDFType == TYPE::CUMULATIVE_VOLUME_DISTRIBUTION)
        convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
    else
        logger(ERROR, "Wrong CDFType");
 
    convertProbabilityDensityToCumulative();
}

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensity(), convertProbabilityDensityToCumulative(), CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_NUMBER_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, ERROR, logger, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, and PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

Referenced by cutoffCumulativeDistribution(), getAreaDx(), getLengthDx(), and getVolumeDx().

convertCumulativeToCumulativeNumberDistribution()

void PSD::convertCumulativeToCumulativeNumberDistribution

(

TYPE

CDFType

)

Convert any CDF to a CUMULATIVE_NUMBER_DISTRIBUTION.

Converts any of the CDFTypes to a the default CUMULATIVE_NUMBER_DISTRIBUTION based on their TYPE.

Parameters

[in]

CDFType

Type of the PDF: CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_AREA_DISTRIBUTION or CUMULATIVE_VOLUME_DISTRIBUTION.

{
    // Ignore the default case.
    if (CDFType == TYPE::CUMULATIVE_NUMBER_DISTRIBUTION)
        return;
 
    convertCumulativeToProbabilityDensity();
 
    if (CDFType == TYPE::CUMULATIVE_LENGTH_DISTRIBUTION)
        convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION);
    else if (CDFType == TYPE::CUMULATIVE_AREA_DISTRIBUTION)
        convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION);
    else if (CDFType == TYPE::CUMULATIVE_VOLUME_DISTRIBUTION)
        convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
    else
        logger(ERROR, "Wrong CDFType");
 
    convertProbabilityDensityToCumulative();
}

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_NUMBER_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, ERROR, logger, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, and PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

Referenced by cutoffCumulativeDistribution().

convertCumulativeToProbabilityDensity()

void PSD::convertCumulativeToProbabilityDensity

(

)

Converts a CDF to a PDF by derivation.

Convert any type of CDF to a PDF. Probabilities are derivated for each radius by substracting the CDF probabilities (i.e. pPDF_i = pCDF_i - pCDF_i-1).

{
    // subtract cumulative probabilities
    Mdouble probabilityOld = 0, probabilityCDF;
    for (auto& it : particleSizeDistribution_)
    {
        probabilityCDF = it.probability;
        it.probability -= probabilityOld;
        probabilityOld = probabilityCDF;
    }
    // check whether probability density distribution is valid
    validateProbabilityDensityDistribution();
}

References particleSizeDistribution_, and validateProbabilityDensityDistribution().

Referenced by convertCumulativeNumberDistributionToCumulative(), convertCumulativeToCumulativeNumberDistribution(), getParticleSizeDistributionByType(), getVolumetricMeanRadius(), insertManuallyByVolume(), and setPSDFromVector().

convertPdfPhiToPdfMeter()

std::vector< std::pair< Mdouble, Mdouble > > PSD::convertPdfPhiToPdfMeter	(	Mdouble	meanInPhi,
		Mdouble	standardDeviationInPhi,
		int	numberOfBins
	)

compute central momenta of the user defined PSD. ‍/ void computeCentralMomenta();

     \brief compute raw momenta of the user defined PSD.
    ‍/

void computeRawMomenta();

/*!

/*!

compute standardised momenta of the user defined PSD. ‍/ void computeStandardisedMomenta();

/*!

get momenta of the user defined PSD. ‍/ std::array<Mdouble, 6> getMomenta() const;

/*!

convert the probabilityDensity of diameter in Phi to the probabilityDensity of diameter in Meter.

generates a normal particle size distribution in Phi units and converts it to the distribution in meters.

Parameters

[in]	meanInPhi	Double representing the mean diameter of the particle size distribution in phi units
[in]	standardDeviationInPhi	Double representing the standard deviation of the particle size distribution (diameter) in phi units
[in]	numberOfBins	Integer determining the number of bins (aka. particle size classes or resolution) of the particle size distribution

                                                                                                                                       {
    Mdouble dMin = meanInPhi - 6 * standardDeviationInPhi;
    Mdouble dMax = meanInPhi + 6 * standardDeviationInPhi;
    std::vector<Mdouble> diameterPhi = helpers::linspace(dMin, dMax, numberOfBins);
    std::vector<Mdouble> pdfPhi;
    for (int i = 0; i < diameterPhi.size(); i++) {
        Mdouble probability = 1 / (standardDeviationInPhi * sqrt(2.0 * constants::pi)) *
                              exp(-0.5 * pow((diameterPhi[i] - meanInPhi) / standardDeviationInPhi, 2));
        pdfPhi.push_back(probability);
    }
    std::vector<Mdouble> diameter;
    std::vector<Mdouble> pdf(pdfPhi.size());
    for(int i = 0; i < diameterPhi.size(); i++){
        Mdouble D = pow(2, -diameterPhi[i])/1000.0;
        diameter.push_back(D);
    }
    pdf[0] = 0;
    for(int i= 0; i < diameterPhi.size()-1; i++){
        Mdouble area = (diameterPhi[i+1]-diameterPhi[i]) * 0.5*(pdfPhi[i+1]+pdfPhi[i]);
        Mdouble dD = fabs(diameter[i+1]-diameter[i]);
        pdf[i+1] = 2.0*area/dD - pdf[i];
    }
    std::vector< std::pair <Mdouble,Mdouble> > vect;
    for (int i=0; i<diameter.size(); i++)
        vect.push_back(std::make_pair(diameter[i],pdf[i]));
    return vect;
}

References D, Eigen::bfloat16_impl::exp(), boost::multiprecision::fabs(), i, helpers::linspace(), constants::pi, Eigen::bfloat16_impl::pow(), and sqrt().

Referenced by setDistributionPhiNormal().

convertProbabilityDensityNumberDistributionToProbabilityDensity()

void PSD::convertProbabilityDensityNumberDistributionToProbabilityDensity

(

TYPE

PDFType

)

Convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to any PDF.

This is a helper function to have full flexibility in converting from a PROBABILITYDENSITY_NUMBER_DISTRIBUTION type to any other PDF type. This can be useful e.g. when wanting to output the PSD with a certain type.

Parameters

PDFType

The PDF type to convert to.

{
    Mdouble sum = 0;
    switch  (PDFType)
    {
        default:
            logger(ERROR, "Wrong PDFType");
            break;
        case TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability *= 0.5 * (it->internalVariable + (it - 1)->internalVariable);
                // old conversion
                //p.probability *= p.internalVariable;
                // sum up probabilities
                sum += it->probability;
            }
            break;
        case TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability *=
                        1.0 / 3.0 * (square(it->internalVariable) + it->internalVariable * (it - 1)->internalVariable +
                                     square((it - 1)->internalVariable));
                // old conversion
                //p.probability *= square(p.internalVariable);
                // sum up probabilities
                sum += it->probability;
            }
            break;
        case TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability *=
                        0.25 * (square(it->internalVariable) + square((it - 1)->internalVariable)) *
                        (it->internalVariable + (it - 1)->internalVariable);
                // old conversion
                //p.probability *= cubic(p.internalVariable);
                // sum up probabilities
                sum += it->probability;
            }
            break;
    }
 
    // normalize
    for (auto& p : particleSizeDistribution_)
    {
        p.probability /= sum;
    }
}

References ERROR, logger, p, particleSizeDistribution_, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, and Eigen::square().

Referenced by convertCumulativeNumberDistributionToCumulative(), convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), getParticleSizeDistributionByType(), getVolumetricMeanRadius(), and insertManuallyByVolume().

convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution()

MERCURYDPM_DEPRECATED void PSD::convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution

(

)

convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

Deprecated:

Use convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION) instead.

converts a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION. Used for getVolumetricMeanRadius() and insertManuallyByVolume().

{
    convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
}

References convertProbabilityDensityNumberDistributionToProbabilityDensity(), and PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

convertProbabilityDensityToCumulative()

void PSD::convertProbabilityDensityToCumulative

(

)

Converts a PDF to a CDF by integration.

Convert any type of PDF to a CDF. Probabilities are integrated for each radius by cumulatively summing them up (i.e. p_i = p_i + p_i-1).

{
    // add up probabilities
    for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
    {
        it->probability += (it - 1)->probability;
    }
    // check whether cumulative distribution is valid
    validateCumulativeDistribution();
}

References particleSizeDistribution_, and validateCumulativeDistribution().

Referenced by convertCumulativeNumberDistributionToCumulative(), convertCumulativeToCumulativeNumberDistribution(), getParticleSizeDistributionByType(), insertManuallyByVolume(), setDistributionPhiNormal(), and setPSDFromVector().

convertProbabilityDensityToProbabilityDensityNumberDistribution()

void PSD::convertProbabilityDensityToProbabilityDensityNumberDistribution

(

TYPE

PDFType

)

Convert any PDF to a PROBABILITYDENSITY_NUMBER_DISTRIBUTION.

Converts any of the PDFTypes to a PROBABILITYDENSITY_NUMBER_DISTRIBUTION based on their TYPE. This is a helper function which enables the convertProbabilityDensityToCumulative function to convert the psd into the default TYPE (CUMULATIVE_NUMBER_DISTRIBUTION).

Parameters

[in]

PDFType

Type of the PDF: PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_AREA_DISTRIBUTION or PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

{
    Mdouble sum = 0;
    switch  (PDFType)
    {
        default:
            logger(ERROR, "Wrong PDFType");
            break;
        case TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability /= 0.5 * (it->internalVariable + (it - 1)->internalVariable);
                // old conversion
                //p.probability /= p.internalVariable;
                // sum up probabilities
                sum += it->probability;
            }
            break;
        case TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability /=
                        1.0 / 3.0 * (square(it->internalVariable) + it->internalVariable * (it - 1)->internalVariable +
                                     square((it - 1)->internalVariable));
                // old conversion
                //p.probability /= square(p.internalVariable);
                // sum up probabilities
                sum += it->probability;
            }
            break;
        case TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability /=
                        0.25 * (square(it->internalVariable) + square((it - 1)->internalVariable)) *
                        (it->internalVariable + (it - 1)->internalVariable);
                // old conversion
                //p.probability /= cubic(p.internalVariable);
                // sum up probabilities
                sum += it->probability;
            }
            break;
    }
 
    // normalize
    for (auto& p : particleSizeDistribution_)
    {
        p.probability /= sum;
    }
}

References ERROR, logger, p, particleSizeDistribution_, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, and Eigen::square().

Referenced by convertCumulativeToCumulativeNumberDistribution(), insertManuallyByVolume(), and setPSDFromVector().

copy()

PSD * PSD::copy

(

)

const

Creates a copy on the heap and returns a pointer.

Copy method; creates a copy on the heap and returns its pointer.

Returns

pointer to the copy on the heap

{
#ifdef DEBUG_CONSTRUCTOR
    std::cout << "PSD::copy() const finished" << std::endl;
#endif
    return new PSD(*this);
}

References PSD().

cutHighSizeRatio()

void PSD::cutHighSizeRatio

(

)

Check if the size ratio is too high and cut it.

Checks if the Size ratio of the PSD is too high and cuts the PSD at head and tail by 10 percent to avoid inaccurate results.

{
    Mdouble SR = getSizeRatio();
    if (SR > 100)
    {
        logger(WARN, "Size ratio > 100; Cutting the PSD to avoid inaccurate results");
        cutoffCumulativeDistribution(TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, 0.1, 0.9, 0.5);
    }
}

References CUMULATIVE_NUMBER_DISTRIBUTION, cutoffCumulativeDistribution(), getSizeRatio(), logger, and WARN.

cutoffAndSqueezeCumulative()

MERCURYDPM_DEPRECATED void PSD::cutoffAndSqueezeCumulative ( Mdouble  quantileMin, Mdouble  quantileMax, Mdouble  squeeze, Mdouble  minPolydispersity = 0.1  )

inline

cutoff the PSD at given quantiles and make it less polydisperse by squeezing it.

Deprecated:

Instead use cutoffAndSqueezeCumulativeDistribution(TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, quantileMin, quantileMax, squeeze, minPolydispersity)

    {
        cutoffAndSqueezeCumulativeDistribution(TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, quantileMin, quantileMax, squeeze, minPolydispersity);
    }

References CUMULATIVE_NUMBER_DISTRIBUTION, and cutoffAndSqueezeCumulativeDistribution().

cutoffAndSqueezeCumulativeDistribution()

void PSD::cutoffAndSqueezeCumulativeDistribution	(	PSD::TYPE	CDFType,
		Mdouble	quantileMin,
		Mdouble	quantileMax,
		Mdouble	squeeze,
		Mdouble	minPolydispersity = 0.1
	)

Cut off the PSD at the given quantiles, with the PSD in the provided CDF form, and make it less polydisperese by squeezing it.

Cuts off the CDF at given minimum and maximum quantiles, applies a minimum polydispersity at the base and squeezes the distribution to make it less polydisperse.

Parameters

[in]	CDFType	The CDF form the PSD should have when applying the cut off and squeeze.
[in]	quantileMin	Undersize quantile to cut off the lower part of the CDF.
[in]	quantileMax	Oversize quantile to cut off the upper part of the CDF.
[in]	squeeze	Applies a squeezing factor ([0,1]) which determines the degree the PDF gets squeezed.
[in]	minPolydispersity	Applies a minimum polydispersity ([0,1]) at the base of the CDF.

{
    Mdouble r50;
    if (CDFType == TYPE::CUMULATIVE_NUMBER_DISTRIBUTION)
        r50 = 0.5 * PSD::getNumberDx(50);
    else if (CDFType == TYPE::CUMULATIVE_LENGTH_DISTRIBUTION)
        r50 = 0.5 * PSD::getLengthDx(50);
    else if (CDFType == TYPE::CUMULATIVE_AREA_DISTRIBUTION)
        r50 = 0.5 * PSD::getAreaDx(50);
    else if (CDFType == TYPE::CUMULATIVE_VOLUME_DISTRIBUTION)
        r50 = 0.5 * PSD::getVolumeDx(50);
    else
        logger(ERROR, "Wrong CDFType");
 
    // cut off
    cutoffCumulativeDistribution(CDFType, quantileMin, quantileMax, minPolydispersity);
    // squeeze psd
    for (auto& p: particleSizeDistribution_)
    {
        p.internalVariable = r50 + (p.internalVariable - r50) * squeeze;
    }
}

References CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_NUMBER_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, cutoffCumulativeDistribution(), ERROR, getAreaDx(), getLengthDx(), getNumberDx(), getVolumeDx(), logger, p, and particleSizeDistribution_.

Referenced by cutoffAndSqueezeCumulative().

cutoffCumulativeDistribution()

void PSD::cutoffCumulativeDistribution	(	PSD::TYPE	CDFType,
		Mdouble	quantileMin,
		Mdouble	quantileMax,
		Mdouble	minPolydispersity = 0.1
	)

Cut off the PSD at the given quantiles, with the PSD in the provided CDF form.

Cuts off the CDF at given minimum and maximum quantiles and applies a minimum polydispersity at the base.

Parameters

[in]	CDFType	The CDF form the PSD should have when applying the cut off.
[in]	quantileMin	Undersize quantile to cut off the lower part of the CDF.
[in]	quantileMax	Oversize quantile to cut off the upper part of the CDF.
[in]	minPolydispersity	Applies a minimum polydispersity ([0,1]) at the base of the CDF.

{
    // Convert the default CUMULATIVE_NUMBER_DISTRIBUTION to the provided CDF type.
    convertCumulativeNumberDistributionToCumulative(CDFType);
 
    Mdouble radiusMin = getRadiusByQuantile(quantileMin);
    Mdouble radiusMax = getRadiusByQuantile(quantileMax);
    // to get a minimum polydispersity at the base
    Mdouble radiusMinCut = std::min(radiusMin * (1 + minPolydispersity), radiusMax);
    // cut off min
    while (particleSizeDistribution_.front().internalVariable <= radiusMinCut)
    {
        particleSizeDistribution_.erase(particleSizeDistribution_.begin());
    }
    particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMinCut, quantileMin});
    particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMin, 0});
    // cut off max
    while (particleSizeDistribution_.back().internalVariable >= radiusMax)
    {
        particleSizeDistribution_.pop_back();
    }
    Mdouble radiusMaxCut = std::max(radiusMax - (1 - minPolydispersity) * (radiusMax - particleSizeDistribution_.back()
            .internalVariable), radiusMin);
    particleSizeDistribution_.push_back({radiusMaxCut, quantileMax});
    particleSizeDistribution_.push_back({radiusMax, 1});
 
    // Convert the provided CDF type back to the default CUMULATIVE_NUMBER_DISTRIBUTION.
    convertCumulativeToCumulativeNumberDistribution(CDFType);
 
    logger(INFO, "PSD was cut successfully and now has a size ratio of %", getSizeRatio());
}

References convertCumulativeNumberDistributionToCumulative(), convertCumulativeToCumulativeNumberDistribution(), getRadiusByQuantile(), getSizeRatio(), INFO, logger, max, min, and particleSizeDistribution_.

Referenced by cutHighSizeRatio(), cutoffAndSqueezeCumulativeDistribution(), and cutoffCumulativeNumber().

cutoffCumulativeNumber()

MERCURYDPM_DEPRECATED void PSD::cutoffCumulativeNumber ( Mdouble  quantileMin, Mdouble  quantileMax, Mdouble  minPolydispersity = 0.1  )

inline

cutoff the PSD at given quantiles.

Deprecated:

Instead use cutoffCumulativeDistribution(PSD::TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, quantileMin, quantileMax, minPolydispersity)

     {
         cutoffCumulativeDistribution(TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, quantileMin, quantileMax, minPolydispersity);
     }

References CUMULATIVE_NUMBER_DISTRIBUTION, and cutoffCumulativeDistribution().

decrementNParticlesPerClass()

void PSD::decrementNParticlesPerClass

(

)

Decrement nParticlesPerClass_ counter.

Decrement the nParticlesPerClass_ counter of the last inserted class by 1. This is utilized when a particle is generated but failed to be placed into an insertionBoundary.

{
    --nParticlesPerClass_[index_];
}

References index_, and nParticlesPerClass_.

decrementVolumePerClass()

void PSD::decrementVolumePerClass

(

Mdouble

volume

)

Decrement volumePerClass_ counter.

Decrement the volumePerClass_ counter of the last inserted class by the volume of the particle. This is utilized when a particle is generated but failed to be placed into an insertionBoundary.

Parameters

[in]

volume

Volume of the particle which failed to be inserted.

{
    volumePerClass_[index_]-= volume;
}

References index_, and volumePerClass_.

drawSample()

Mdouble PSD::drawSample

(

)

Draw a sample radius from a CUMULATIVE_NUMBER_DISTRIBUTION.

Draws a sample probability of a real uniform distribution. A random number is generated in range [0,1] and by linear interpolation a suitable radius is returned. This function is only valid for CumulativeNumberDistributions as particles are drawn and not volumes, areas or lengths.

Returns

A Radius for a randomly assigned probability.

{
    // draw a number between 0 and 1, uniformly distributed
    Mdouble prob = random_.getRandomNumber(0, 1);
    // find the interval [low,high] of the psd in which the number sits
    auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), prob);
    auto low = std::max(particleSizeDistribution_.begin(), high - 1);
    Mdouble rMin = low->internalVariable;
    Mdouble rMax = high->internalVariable;
    Mdouble pMin = low->probability;
    Mdouble pMax = high->probability;
    // interpolate linearly between [low.radius,high.radius] (assumption: CDF is piecewise linear)
    Mdouble a = (prob - pMin) / (pMax - pMin);
    return a * rMin + (1 - a) * rMax;
}

References a, RNG::getRandomNumber(), max, particleSizeDistribution_, and random_.

getAreaDx()

Mdouble PSD::getAreaDx

(

Mdouble

x

)

const

Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the area based PSD.

Gets the diameter from a certain percentile of the area based PSD.

Returns

A double which is the diameter corresponding to a certain percentile.

Parameters

[in]

x

double [0, 100] which determines the obtained diameter as a percentile of the PSD.

{
    logger.assert_always(x >= 0 && x <= 100, "Percentile % is not between 0 and 100.", x);
    PSD psd = *this;
    psd.convertCumulativeNumberDistributionToCumulative(TYPE::CUMULATIVE_AREA_DISTRIBUTION);
    return 2.0 * psd.getRadiusByQuantile(x / 100);
}

References convertCumulativeNumberDistributionToCumulative(), CUMULATIVE_AREA_DISTRIBUTION, getRadiusByQuantile(), logger, and plotDoE::x.

Referenced by cutoffAndSqueezeCumulativeDistribution().

getDistributionLogNormal()

static PSD PSD::getDistributionLogNormal ( Mdouble  mean, Mdouble  standardDeviation, int  numberOfBins  )

inlinestatic

                                                                                                   {
        PSD psd;
        psd.setDistributionLogNormal(mean, standardDeviation, numberOfBins);
        return psd;
    }

References setDistributionLogNormal().

Referenced by DistributionSelfTest::setupInitialConditions().

getDistributionNormal()

static PSD PSD::getDistributionNormal ( Mdouble  mean, Mdouble  standardDeviation, int  numberOfBins  )

inlinestatic

                                                                                                {
        PSD psd;
        psd.setDistributionNormal(mean, standardDeviation, numberOfBins);
        return psd;
    }

References setDistributionNormal().

Referenced by Material::setPSD(), Drum::setupInitialConditions(), and RotatingDrumWet::setupInitialConditions().

getInsertedParticleNumber()

int PSD::getInsertedParticleNumber

(

)

const

Get the number of particles already inserted into the simulation.

Gets the number of particles already inserted into the simulation by summing up the particles inserted in each class.

Returns

An integer containing the number of particles already inserted into the simulation.

{
    int sum = 0;
    for (auto& it: nParticlesPerClass_)
        sum += it;
    return sum;
}

References nParticlesPerClass_.

getLengthDx()

Mdouble PSD::getLengthDx

(

Mdouble

x

)

const

Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the length based PSD.

Gets the diameter from a certain percentile of the length based PSD.

Returns

A double which is the diameter corresponding to a certain percentile.

Parameters

[in]

x

double [0, 100] which determines the obtained diameter as a percentile of the PSD.

{
    logger.assert_always(x >= 0 && x <= 100, "Percentile % is not between 0 and 100.", x);
    PSD psd = *this;
    psd.convertCumulativeNumberDistributionToCumulative(TYPE::CUMULATIVE_LENGTH_DISTRIBUTION);
    return 2.0 * psd.getRadiusByQuantile(x / 100);
}

References convertCumulativeNumberDistributionToCumulative(), CUMULATIVE_LENGTH_DISTRIBUTION, getRadiusByQuantile(), logger, and plotDoE::x.

Referenced by cutoffAndSqueezeCumulativeDistribution().

getMaxRadius()

Mdouble PSD::getMaxRadius

(

)

const

Get largest radius of the PSD.

Gets the maximum radius of the PSD.

Returns

A double which corresponds to the maximum radius of the PSD.

{
    return particleSizeDistribution_.back().internalVariable;
}

References particleSizeDistribution_.

Referenced by getSizeRatio(), GranuDrum::GranuDrum(), and GranuHeap::GranuHeap().

getMinRadius()

Mdouble PSD::getMinRadius

(

)

const

Get smallest radius of the PSD.

Gets the minimal radius of the PSD.

Returns

A double which corresponds to the minimal radius of the PSD.

{
    return particleSizeDistribution_[0].internalVariable;
}

References particleSizeDistribution_.

Referenced by getSizeRatio(), GranuHeap::GranuHeap(), Calibration::setSpecies(), and Material::setSpecies().

getNumberDx()

Mdouble PSD::getNumberDx

(

Mdouble

x

)

const

Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the number based PSD.

Gets the diameter from a certain percentile of the number based PSD.

Returns

A double which is the diameter corresponding to a certain percentile.

Parameters

[in]

x

double [0, 100] which determines the obtained diameter as a percentile of the PSD.

{
    logger.assert_always(x >= 0 && x <= 100, "Percentile % is not between 0 and 100.", x);
    return 2.0 * getRadiusByQuantile(x / 100);
}

References getRadiusByQuantile(), logger, and plotDoE::x.

Referenced by cutoffAndSqueezeCumulativeDistribution(), and GranuHeap::GranuHeap().

getParticleSizeDistribution()

std::vector< DistributionElements > PSD::getParticleSizeDistribution

(

)

const

Get the PSD vector.

Gets the vector containing radii and probabilities of the PSD.

Returns

A vector containing radii and probabilities of the PSD.

{
    return particleSizeDistribution_;
}

References particleSizeDistribution_.

Referenced by getParticleSizeDistributionByType(), DistributionPhiNormalSelfTest::setupInitialConditions(), and SetDistributionPhiNormalSelfTest::setupInitialConditions().

getParticleSizeDistributionByType()

std::vector< DistributionElements > PSD::getParticleSizeDistributionByType	(	TYPE	PSDType,
		Mdouble	scalingFactor = 1.0
	)		const

Gets the PSD vector, converted to the preferred type.

This returns the PSD vector converted to the specified type, with an optional scale factor for the internal variable. This can be used e.g. to output the PSD with a certain type.

Parameters

PSDType	The type the PSD should be.
scalingFactor	The scale factor by which the the internal variable is multiplied.

Returns

The converted PSD vector.

{
    PSD psd = *this;
 
    switch (PSDType)
    {
        // default case is a CUMULATIVE_NUMBER_DISTRIBUTION
        default:
            break;
        case TYPE::CUMULATIVE_LENGTH_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION);
            psd.convertProbabilityDensityToCumulative();
            break;
        case TYPE::CUMULATIVE_AREA_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION);
            psd.convertProbabilityDensityToCumulative();
            break;
        case TYPE::CUMULATIVE_VOLUME_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
            psd.convertProbabilityDensityToCumulative();
            break;
        case TYPE::PROBABILITYDENSITY_NUMBER_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            break;
        case TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(PSDType);
            break;
        case TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(PSDType);
            break;
        case TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            psd.convertCumulativeToProbabilityDensity();
            psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(PSDType);
            break;
    }
 
    auto psdVector = psd.getParticleSizeDistribution();
    for (auto& de : psdVector)
        de.internalVariable *= scalingFactor;
 
    return psdVector;
}

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensity(), convertProbabilityDensityToCumulative(), CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, getParticleSizeDistribution(), PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_NUMBER_DISTRIBUTION, and PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

Referenced by ParticleHandler::saveNumberPSDtoCSV(), and testPSD().

getQuantileByRadius()

Mdouble PSD::getQuantileByRadius

(

Mdouble

radius

)

const

Calculates the quantile corresponding to a certain radius.

Calculates the quantile corresponding to a certain radius.

Parameters

radius

The radius to find the quantile for.

Returns

The quantile found.

{
    auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(),
                                 radius, [](const DistributionElements& rp, Mdouble value){ return rp.internalVariable < value; });
    auto low = std::max(particleSizeDistribution_.begin(), high - 1);
    if (high->internalVariable == low->internalVariable)
        return high->probability;
    else
        return low->probability + (high->probability - low->probability) * (radius - low->internalVariable) / (high->internalVariable - low->internalVariable);
}

References max, particleSizeDistribution_, UniformPSDSelfTest::radius, and Eigen::value.

getRadius()

Mdouble PSD::getRadius

(

int

index

)

const

get a radius at a certain index of the particleSizeDistribution_

 \details Compute the raw momenta of the inserted PSD by converting it to a PDF, calculating the moments m_i according
 to \f$ m_i = \sum_{j=0}^n \f$ scaling it by the number
 of particles inserted into the simulation (moments are computed from a PROBABILITYDENSITY_NUMBER_DISTRIBUTION) and
 converting it back to a CDF.
 See <a href="https://en.wikipedia.org/wiki/Moment_(mathematics)#Significance_of_the_moments">Wikipedia</a> for
 details.
&zwj;/

void PSD::computeRawMomenta() { convertCumulativeToProbabilityDensity(); for (size_t im = 0; im < momenta_.size(); ++im) { // prevent summing up of moments for each time step of the simulation momenta_[im] = 0; for (auto& it : particleSizeDistribution_) { momenta_[im] += std::pow(it.internalVariable, im) * it.probability; } } // zeroth-moment equals particle number for a PROBABILITYDENSITY_NUMBER_DISTRIBUTION momenta_[0] *= getInsertedParticleNumber(); convertProbabilityDensityToCumulative(); }

/*!

Compute the central momenta of the PSD from their respective raw momenta. ‍/ void PSD::computeCentralMomenta() { computeRawMomenta(); Mdouble mean = momenta_[1]; momenta_[5] += -5 * mean * momenta_[4] + 10 * mean * mean * momenta_[3]

• 10 * mean * mean * mean * momenta_[2] + 4 * mean * mean * mean * mean * mean; momenta_[4] += -4 * mean * momenta_[3] + 6 * mean * mean * momenta_[2] - 3 * mean * mean * mean * mean; momenta_[3] += -3 * mean * momenta_[2] + 2 * mean * mean * mean; momenta_[2] += -mean * mean; }

/*!

Compute the standardised momenta of the PSD from their respective central momenta. ‍/ void PSD::computeStandardisedMomenta() { computeCentralMomenta(); Mdouble std = std::sqrt(momenta_[2]); momenta_[3] /= std * std * std; momenta_[4] /= std * std * std * std; momenta_[5] /= std * std * std * std * std; }

/*!

Gets the radius of a particle from the PSD.

Parameters

[in]

i

An integer which is the index of the vector containing the radius.

Returns

A double which corresponds to the radius of the particle.

{
    return particleSizeDistribution_[index].internalVariable;
}

References particleSizeDistribution_.

getRadiusByQuantile()

Mdouble PSD::getRadiusByQuantile

(

Mdouble

quantile

)

const

Calculate the quantile of the PSD.

gets the radius from a certain quantile of the PSD

Returns

A double which is the radius corresponding to a certain quantile of the PSD.

Parameters

[in]

quantile

double which determines the returned radius as a quantile of the PSD.

{
    logger.assert_always(quantile <= 1 && quantile >= 0, "quantile is not between 0 and 1");
    // find the quantile corresponding to the PSD
    auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), quantile);
    auto low = std::max(particleSizeDistribution_.begin(), high - 1);
    if (high->probability == low->probability)
        return high->internalVariable;
    else
        return low->internalVariable +
               (high->internalVariable - low->internalVariable) * (quantile - low->probability) /
               (high->probability - low->probability);
}

References logger, max, and particleSizeDistribution_.

Referenced by cutoffCumulativeDistribution(), getAreaDx(), getLengthDx(), getNumberDx(), and getVolumeDx().

getSizeRatio()

Mdouble PSD::getSizeRatio

(

)

const

get the size ratio (width) of the PSD.

Gets the size ratio (width) of the PSD defined by the ratio of minimum to maximum particle radius.

Returns

A double which corresponds to the size ratio of the PSD.

{
    Mdouble rMin = getMaxRadius();
    Mdouble rMax = getMinRadius();
    return rMin / rMax;
}

References getMaxRadius(), and getMinRadius().

Referenced by cutHighSizeRatio(), cutoffCumulativeDistribution(), and setPSDFromCSV().

getVolumeDx()

Mdouble PSD::getVolumeDx

(

Mdouble

x

)

const

Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the volume based PSD.

Gets the diameter from a certain percentile of the volume based PSD.

Returns

A double which is the diameter corresponding to a certain percentile.

Parameters

[in]

x

double [0, 100] which determines the obtained diameter as a percentile of the PSD.

{
    logger.assert_always(x >= 0 && x <= 100, "Percentile % is not between 0 and 100.", x);
    PSD psd = *this;
    psd.convertCumulativeNumberDistributionToCumulative(TYPE::CUMULATIVE_VOLUME_DISTRIBUTION);
    return 2.0 * psd.getRadiusByQuantile(x / 100);
}

References convertCumulativeNumberDistributionToCumulative(), CUMULATIVE_VOLUME_DISTRIBUTION, getRadiusByQuantile(), logger, and plotDoE::x.

Referenced by cutoffAndSqueezeCumulativeDistribution(), GranuDrum::GranuDrum(), GranuHeap::GranuHeap(), Calibration::setSpecies(), and Material::setSpecies().

getVolumetricMeanRadius()

Mdouble PSD::getVolumetricMeanRadius

(

)

const

get a volumetric mean radius of the PSD.

Gets a radius such that a monodisperse system has the same number of particles as a polydisperse system. (i.e. mean += p_i * 0.5*(r_i^3 + r_i-1^3)

Returns

A double which corresponds to the volumetric mean radius.

{
    PSD psd = *this;
    psd.convertCumulativeToProbabilityDensity();
    psd.convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
    Mdouble mean = 0;
    for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
        mean += it->probability * 0.5 * (cubic(it->internalVariable) + cubic((it - 1)->internalVariable));
    mean = pow(mean, 1. / 3.);
    return mean;
}

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensity(), mathsFunc::cubic(), particleSizeDistribution_, Eigen::bfloat16_impl::pow(), and PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

insertManuallyByVolume()

Mdouble PSD::insertManuallyByVolume

(

Mdouble

volume

)

Draw sample radius manually per size class and check the volumeAllowed of each size class to insert the PSD as accurate as possible.

Draws a sample probability of a real uniform distribution within a given size class. A random number is generated in the size class range [r_i,r_i+1] and by linear interpolation a suitable radius is returned. from a CUMULATIVE_VOLUME_DISTRIBUTION the volumeAllowed of each class is checked to insert the PSD as accurate as possible. This function is only valid for cumulativeNumberDistributions as particles are drawn and not volumes, areas or lengths. Furthermore this insertion routine is most accurate for non-continuous particle insertion.

Parameters

[in]

volume

volume of the geometry to be filled.

Returns

A Radius for a randomly assigned probability of a specific size class.

{
    // initialize the particleNumberPerClass vector if empty.
    if (nParticlesPerClass_.empty())
        nParticlesPerClass_.resize(particleSizeDistribution_.size());
    // initialize the particleNumberPerClass vector if empty.
    if (volumePerClass_.empty())
        volumePerClass_.resize(particleSizeDistribution_.size());
 
    for (auto it = particleSizeDistribution_.end() - 1; it != particleSizeDistribution_.begin(); --it)
    {
        static std::mt19937 gen(0);
        static std::uniform_real_distribution<Mdouble> dist(0, 1);
        Mdouble prob = dist(gen);
        // map the probability to the current class interval
        prob = (it - 1)->probability + (it->probability - (it - 1)->probability) * prob;
        // find the interval [low,high] of the psd in which the number sits
        auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), prob);
        auto low = std::max(particleSizeDistribution_.begin(), high - 1);
        Mdouble rMin = low->internalVariable;
        Mdouble rMax = high->internalVariable;
        Mdouble pMin = low->probability;
        Mdouble pMax = high->probability;
        // interpolate linearly between [low.internalVariable,high.internalVariable] (assumption: CDF is piecewise linear)
        index_ = std::distance(particleSizeDistribution_.begin(), high);
        Mdouble a = (prob - pMin) / (pMax - pMin);
        Mdouble rad = a * rMin + (1 - a) * rMax;
        //\todo generalize this to work for non-spherical particles
        // Compare inserted volume against volumeAllowed (only works for spherical at the moment)
        Mdouble radVolume = 4.0 / 3.0 * constants::pi * rad * rad * rad;
        convertCumulativeToProbabilityDensity();
        convertProbabilityDensityNumberDistributionToProbabilityDensity(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
        Mdouble volumeAllowed = particleSizeDistribution_[index_].probability * volume;
        convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
        convertProbabilityDensityToCumulative();
        if (volumePerClass_[index_] > volumeAllowed)
        {
            continue;
        }
        else if (radVolume + volumePerClass_[index_] > volumeAllowed)
        {
            Mdouble differenceVolumeLow = -(volumePerClass_[index_] - volumeAllowed);
            Mdouble differenceVolumeHigh = radVolume + volumePerClass_[index_] - volumeAllowed;
            Mdouble volumeRatio = differenceVolumeLow / differenceVolumeHigh;
            // If volumeRatio > 1 it will insert anyways, because it should be no problem for the distribution
            prob = dist(gen);
            if (prob <= volumeRatio)
            {
                ++nParticlesPerClass_[index_];
                volumePerClass_[index_] += radVolume;
                return rad;
            }
            else
            {
                continue;
            }
        }
        else
        {
            ++nParticlesPerClass_[index_];
            volumePerClass_[index_] += radVolume;
            return rad;
        }
        
    }
    return 0;
}

References a, convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensity(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), index_, max, nParticlesPerClass_, particleSizeDistribution_, constants::pi, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, and volumePerClass_.

linspace()

static std::vector<Mdouble> PSD::linspace ( Mdouble  Min, Mdouble  Max, int  numberOfBins  )

static

create a vector of linearly spaced values.

printPSD()

void PSD::printPSD

(

)

const

Prints radii and probabilities of the PSD vector.

Prints the radii [m] and probabilities [%] of the psd vector. It currently supports nanometers, micrometers, milimeters and meters as size units.

{
    if (particleSizeDistribution_.front().internalVariable > 1e-1)
    {
        for (const auto p: particleSizeDistribution_)
        {
            logger(INFO, "%m\t %\%", p.internalVariable, p.probability * 100);
        }
    }
    else if (particleSizeDistribution_.front().internalVariable > 1e-4)
    {
        for (const auto p: particleSizeDistribution_)
        {
            logger(INFO, "%mm\t %\%", p.internalVariable * 1000, p.probability * 100);
        }
    }
    else if (particleSizeDistribution_.front().internalVariable > 1e-7)
    {
        for (const auto p: particleSizeDistribution_)
        {
            logger(INFO, "%um\t %\%", p.internalVariable * 1e6, p.probability * 100);
        }
    }
    else
    {
        for (const auto p: particleSizeDistribution_)
        {
            logger(INFO, "%nm\t %\%", p.internalVariable * 1e9, p.probability * 100);
        }
    }
}

References e(), INFO, logger, p, and particleSizeDistribution_.

scaleParticleSize()

void PSD::scaleParticleSize

(

double

scale

)

Scales all particle sizes by a factor.

{
    for (auto& p : particleSizeDistribution_) {
        p.internalVariable *= scale;
    }
}

References p, and particleSizeDistribution_.

Referenced by scaleParticleSizeAuto().

scaleParticleSizeAuto()

double PSD::scaleParticleSizeAuto	(	int	numberOfParticles,
		double	targetVolume,
		bool	allowScaleDown = false
	)

Scales all particle sizes, such that the total volume of N particles approximately equals the target volume.

Parameters

numberOfParticles	The number of particles N.
targetVolume	The volume to match.
allowScaleDown	Whether or not to scale down the particle sizes when the scale factor is smaller than 1.

Returns

The scale factor used.

{
    // Calculate the mean particle volume, using the Mean Value Theorem for Integrals. I.e. the mean particle volume
    // over an infinite number of particles.
    double meanParticleVolume = 0.0;
    for (auto it = particleSizeDistribution_.begin(); it < particleSizeDistribution_.end() - 1; it++)
    {
        double radiusLeft = it->internalVariable;
        double radiusRight = (it+1)->internalVariable;
        double probabilityLeft = it->probability;
        double probabilityRight = (it+1)->probability;
 
        double meanVolume;
        if (radiusLeft != radiusRight)
        {
            // The CDF is assumed piecewise linear, so the radius is linearly interpolated.
            // Integrate the volume from radius left to right, i.e. integral of 4/3 pi r^3 gives 1/3 pi (r_r^4 - r_l^4).
            double volumeIntegral = (std::pow(radiusRight, 4) - std::pow(radiusLeft, 4)) * constants::pi / 3.0;
            // Calculate the mean volume that represents this bin.
            meanVolume = volumeIntegral / (radiusRight - radiusLeft);
        }
        else
        {
            meanVolume = std::pow(radiusLeft, 3) * constants::pi * 4.0 / 3.0;
        }
 
        // Calculate the volume that this bin accounts for, and add to total.
        meanParticleVolume += meanVolume * (probabilityRight - probabilityLeft);
    }
 
    // The volume that the number of particles currently fills, and the scale factor to make it match the target volume.
    double totalParticleVolume = numberOfParticles * meanParticleVolume;
    double scale = std::cbrt(targetVolume / totalParticleVolume);
 
    // The number of (unscaled) particles needed to fill the target volume.
    unsigned long long numberOfParticlesExpected = std::pow(scale, 3) * numberOfParticles;
    logger(INFO, "Rescaling the particle size based on number of particles %, and target volume %.", numberOfParticles, targetVolume);
    logger(INFO, "Expected number of particles %, scale factor %.", numberOfParticlesExpected, scale);
 
    // Always scale up, only scale down when requested.
    if (scale > 1.0 || allowScaleDown)
    {
        scaleParticleSize(scale);
        logger(INFO, "Rescaling applied.");
    }
    else
    {
        logger(INFO, "Rescaling not applied, since scale factor is less than or equal to 1.");
        scale = 1.0; // Reassign so that the returned value is correct.
    }
 
    return scale;
}

References Eigen::numext::cbrt(), INFO, logger, particleSizeDistribution_, constants::pi, Eigen::bfloat16_impl::pow(), and scaleParticleSize().

setDistributionLogNormal()

void PSD::setDistributionLogNormal	(	Mdouble	mean,
		Mdouble	standardDeviation,
		int	numberOfBins
	)

create a PSD vector for a normal distribution.

sets the particle size distribution to a discretised lognormal cumulative number distribution, which generates random values of radius following a lognormal distribution and iterates until the mass-based D50 is close to two times the user-defined mean radius.

Parameters

[in]	mean	Double representing the mean of the particle size distribution in meters (radius)
[in]	standardDeviation	Double representing the standard deviation of the particle size distribution in meters (radius)
[in]	numberOfBins	Integer determining the number of bins (aka. particle size classes or resolution) of the particle size distribution See https://en.cppreference.com/w/cpp/numeric/random/lognormal_distribution

{
    //setDistributionNormal(mean, standardDeviation, numberOfBins);
    if (!particleSizeDistribution_.empty())
    {
        particleSizeDistribution_.clear();
    }
    Mdouble radMin = mean - 3 * standardDeviation;
    Mdouble radMax = mean + 3 * standardDeviation;
    std::vector<Mdouble> radii = helpers::linspace(radMin, radMax, numberOfBins);
    std::vector<Mdouble> probabilities;
    for (int i = 0; i < radii.size(); i++)
    {
        Mdouble probability = 0.5 * (1 + erf((radii[i] - mean) / (sqrt(2) * standardDeviation)));
        probabilities.push_back(probability);
    }
    for (int j = 0; j < radii.size(); j++)
    {
        particleSizeDistribution_.push_back({radii[j], probabilities[j]});
    }
    validateCumulativeDistribution();
    for (auto& p: particleSizeDistribution_)
    {
        p.internalVariable = exp(p.internalVariable);
    }
}

References Eigen::bfloat16_impl::exp(), i, j, helpers::linspace(), p, particleSizeDistribution_, sqrt(), and validateCumulativeDistribution().

Referenced by getDistributionLogNormal().

setDistributionNormal()

void PSD::setDistributionNormal	(	Mdouble	mean,
		Mdouble	standardDeviation,
		int	numberOfBins
	)

create a PSD vector for a normal distribution.

sets the particle size distribution to a discretised normal (gaussian) cumulative number distribution, which covers the range of 3 * standardDeviation (99,73% of all values covered).

Parameters

[in]	mean	Double representing the mean of the particle size distribution
[in]	standardDeviation	Double representing the standard deviation of the particle size distribution
[in]	numberOfBins	Integer determining the number of bins (aka. particle size classes or resolution) of the particle size distribution

{
    logger.assert_always(mean > 3 * standardDeviation,
                         "Reduce standardDeviation of your normal distribution to avoid negative radii; The mean "
                         "should be greater than 3*standardDeviation");
    if (!particleSizeDistribution_.empty())
    {
        particleSizeDistribution_.clear();
    }
    Mdouble radMin = mean - 3 * standardDeviation;
    Mdouble radMax = mean + 3 * standardDeviation;
    std::vector<Mdouble> radii = helpers::linspace(radMin, radMax, numberOfBins);
    std::vector<Mdouble> probabilities;
    for (int i = 0; i < radii.size(); i++)
    {
        Mdouble probability = 0.5 * (1 + erf((radii[i] - mean) / (sqrt(2) * standardDeviation)));
        probabilities.push_back(probability);
    }
    for (int j = 0; j < radii.size(); j++)
    {
        particleSizeDistribution_.push_back({radii[j], probabilities[j]});
    }
    validateCumulativeDistribution();
}

References i, j, helpers::linspace(), logger, particleSizeDistribution_, sqrt(), and validateCumulativeDistribution().

Referenced by getDistributionNormal(), and DistributionToPSDSelfTest::setupInitialConditions().

setDistributionPhiNormal()

void PSD::setDistributionPhiNormal	(	Mdouble	D50,
		Mdouble	standardDeviationInPhi,
		int	numberOfBins
	)

create a PSD vector for a normal distribution in Phi Units that has the demanded D50.

sets the particle size distribution to a discretised cumulative number distribution, which is a normal distribution in phi units and has the demanded D50. D50 is the mass-based median grain size, while the Phi scale is a logarithmic scale computed by: Phi = -log2(1000*D). It can be used for visualizing sediment grain size distributions over a wide range of grain sizes. The combination of D50 and standardDeviation in Phi units is a representative measure of grain size distribution in the field of Sedimentology. Source: Donoghue, J.F. (2016). Phi Scale. In: Kennish, M.J. (eds) Encyclopedia of Estuaries. Encyclopedia of Earth Sciences Series. Springer, Dordrecht. https://doi.org/10.1007/978-94-017-8801-4_277

Parameters

[in]	D50	Double representing the mass-based D50 of the particle size distribution in meters
[in]	standardDeviationInPhi	Double representing the standard deviation of the particle size distribution (diameter) in phi units
[in]	numberOfBins	Integer determining the number of bins (aka. particle size classes or resolution) of the particle size distribution

                                                                                               {
    Mdouble PhiMean = -log2(1000*D50);
    std::vector<std::pair<Mdouble,Mdouble>> vect = convertPdfPhiToPdfMeter(PhiMean, standardDeviationInPhi, numberOfBins);
    std::sort(vect.begin(), vect.end());
    Mdouble D50Test = computeD50(vect);
    Mdouble res = D50Test - D50;
    Mdouble tol = D50*1e-5;
    int itter = 0;
    Mdouble dPhiMean, slope, res_old;
    while (res > tol && itter < 100) {
        if(itter==0) dPhiMean = 0.1;
        else {
            slope = (res - res_old)/dPhiMean;
            dPhiMean = -res/slope;
        }
        PhiMean = PhiMean+dPhiMean;
        vect = convertPdfPhiToPdfMeter(PhiMean, standardDeviationInPhi, numberOfBins);
        std::sort(vect.begin(), vect.end());
        D50Test = computeD50(vect);
        res_old = res;
        res = D50Test - D50;
        itter++;
    }
    logger(INFO, "D50 %, Final D50Test %", D50, D50Test);
 
    Mdouble probTot=0;
    for(int i = 0; i < vect.size(); i++){
        probTot += vect[i].second;
    }
    for (int i = 0; i < vect.size(); i++)
    {
        vect[i].second = vect[i].second/probTot;
    }
 
    if (!particleSizeDistribution_.empty())
    {
        particleSizeDistribution_.clear();
    }
    for (int j = 0; j < vect.size(); j++)
    {
        particleSizeDistribution_.push_back({vect[j].first*0.5, vect[j].second});
    }
    convertProbabilityDensityToCumulative();
    validateCumulativeDistribution();
}

References computeD50(), convertPdfPhiToPdfMeter(), convertProbabilityDensityToCumulative(), e(), i, INFO, j, log2(), logger, particleSizeDistribution_, res, and validateCumulativeDistribution().

Referenced by DistributionPhiNormalSelfTest::setupInitialConditions(), and SetDistributionPhiNormalSelfTest::setupInitialConditions().

setDistributionUniform()

void PSD::setDistributionUniform	(	Mdouble	radMin,
		Mdouble	radMax,
		int	numberOfBins
	)

create a PSD vector for a uniform distribution.

sets the particle size distribution to a discretised uniform (linear) cumulative number distribution

Parameters

[in]	radMin	Double representing The smallest particle radius of the particle size distribution
[in]	radMax	Double representing the biggest particle radius of the particle size distribution
[in]	numberOfBins	Integer determining the number of bins (aka. particle size classes or resolution) of the particle size distribution

{
    if (!particleSizeDistribution_.empty())
    {
        particleSizeDistribution_.clear();
    }
    std::vector<Mdouble> probabilities = helpers::linspace(0.0, 1.0, numberOfBins);
    std::vector<Mdouble> radii = helpers::linspace(radMin, radMax, numberOfBins);
    // combine radii and probabilities
    for (int i = 0; i < radii.size(); ++i)
    {
        particleSizeDistribution_.push_back({radii[i], probabilities[i]});
    }
    validateCumulativeDistribution();
}

References i, helpers::linspace(), particleSizeDistribution_, and validateCumulativeDistribution().

Referenced by BoundariesSelfTest::BoundariesSelfTest(), FluxBoundarySelfTest::FluxBoundarySelfTest(), RandomClusterInsertionBoundary::set(), HopperInsertionBoundary::set(), CylinderInsertionBoundary::set(), SphereInsertionBoundary::set(), ChuteInsertionBoundary::set(), CubeInsertionBoundary::set(), InsertionBoundarySelfTest::setupInitialConditions(), FullRestartTest::setupInitialConditions(), Chute::setupInitialConditions(), ChuteWithHopper::setupInitialConditions(), and T_protectiveWall::T_protectiveWall().

setFixedSeed()

void PSD::setFixedSeed ( int  seed )

inline

set a fixed seed for the random number generator; this is used for the PSDSelfTest to reproduce results

    {
//        random_.setLinearCongruentialGeneratorParmeters(0,0,1);
        random_.setRandomSeed(seed);
    }

References random_, and RNG::setRandomSeed().

setParticleSizeDistribution()

void PSD::setParticleSizeDistribution

(

std::vector< DistributionElements >

particleSizeDistribution

)

set the PSD by a suitable vector.

Deprecated:

Use setPSDFromVector() instead.

sets the PSD from a vector of the DistributionElements class; used to read in PSDs from a restart file.

Parameters

[in]

particleSizeDistribution

vector containing the radii and probabilities specifying the PSD.

{
    particleSizeDistribution_ = particleSizeDistribution;
}

References particleSizeDistribution_.

setPSDFromCSV()

void PSD::setPSDFromCSV	(	const std::string &	fileName,
		TYPE	PSDType,
		bool	headings = false,
		Mdouble	unitScalingFactorRadii = 1.0
	)

read in the PSD vector with probabilities and radii saved in a .csv file.

creates the PSD vector from probabilities and radii saved in a .csv file. Radii should be located at the first column and probabilities at the second column of the file. The Type of PSD will be converted to the default cumulative number distribution function (CUMULATIVE_NUMBER_DISTRIBUTION) for further processing.

Parameters

[in]	fileName	Name of the .csv file containing the radii and probabilities of the PSD.
[in]	PSDType	Type of the PSD: CUMULATIVE_VOLUME_DISTRIBUTION, CUMULATIVE_NUMBER_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_AREA_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, PROBABILITYDENSITY_NUMBER_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION PROBABILITYDENSITY_AREA_DISTRIBUTION.
[in]	headings	If TRUE the file is assumed to have headings and the first row will be skipped. If FALSE the file has no headings and the file will be read in as is. Default is FALSE.
[in]	unitScalingFactorRadii	Scaling factor of radii to match SI-units.

{
    // read in csv file using the csvReader class
    csvReader csv;
    csv.setHeader(headings);
    csv.read(fileName);
    std::vector<Mdouble> radii = csv.getFirstColumn(unitScalingFactorRadii);
    std::vector<Mdouble> probabilities = csv.getSecondColumn(1.0);
    logger.assert_always(radii.size() == probabilities.size(), "The radii and probabilities vector have to be the "
                                                               "same size");
    // combine radii and probabilities
    std::vector<DistributionElements> psd;
    for (int i = 0; i < radii.size(); ++i)
    {
        psd.push_back({radii[i], probabilities[i]});
    }
    logger.assert_always(!psd.empty(), "PSD cannot be empty");
    setPSDFromVector(psd, PSDType);
    logger(INFO, "A PSD with size ratio % is now ready to be set", getSizeRatio());
}

References UniformPSDSelfTest::fileName, csvReader::getFirstColumn(), csvReader::getSecondColumn(), getSizeRatio(), i, INFO, logger, csvReader::read(), csvReader::setHeader(), and setPSDFromVector().

Referenced by MultiplePSDSelfTest::setupInitialConditions(), PSDManualInsertionSelfTest::setupInitialConditions(), and PSDSelfTest::setupInitialConditions().

setPSDFromVector()

void PSD::setPSDFromVector	(	std::vector< DistributionElements >	psdVector,
		TYPE	PSDType
	)

Sets the PSD from a vector with DistributionElements.

creates the psd vector from radii and probabilities filled in by hand. The Type of PSD will be converted to the default cumulative number distribution function (CUMULATIVE_NUMBER_DISTRIBUTION) for further processing.

Parameters

[in]	psdVector	Vector containing radii and probabilities ([0,1]).
[in]	PSDType	Type of the PSD: CUMULATIVE_VOLUME_DISTRIBUTION, CUMULATIVE_NUMBER_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_AREA_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, PROBABILITYDENSITY_NUMBER_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_AREA_DISTRIBUTION.

Todo:

JWB: Is a PSD ever allowed to be empty? drawSample() will cause a segfault in that case. I'm not throwing an error here, as that causes the UnitTests/FullRestartSelfTest to fail, due to the InsertionBoundary being allowed to use (and defaults to using) an empty PSD.

{
    particleSizeDistribution_ = psdVector;
 
    // To prevent errors when the PSD is empty, throw a warning and return.
    if (particleSizeDistribution_.empty())
    {
        logger(WARN, "Setting an empty PSD!");
        return;
    }
 
    switch (PSDType)
    {
        // default case is a CUMULATIVE_NUMBER_DISTRIBUTION
        default:
            validateCumulativeDistribution();
            break;
        case TYPE::CUMULATIVE_LENGTH_DISTRIBUTION:
            validateCumulativeDistribution();
            convertCumulativeToProbabilityDensity();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(
                    TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION);
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::CUMULATIVE_AREA_DISTRIBUTION:
            validateCumulativeDistribution();
            convertCumulativeToProbabilityDensity();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION);
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::CUMULATIVE_VOLUME_DISTRIBUTION:
            validateCumulativeDistribution();
            convertCumulativeToProbabilityDensity();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(
                    TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::PROBABILITYDENSITY_NUMBER_DISTRIBUTION:
            validateProbabilityDensityDistribution();
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION:
            validateProbabilityDensityDistribution();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(PSDType);
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION:
            validateProbabilityDensityDistribution();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(PSDType);
            convertProbabilityDensityToCumulative();
            break;
        case TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            validateProbabilityDensityDistribution();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(PSDType);
            convertProbabilityDensityToCumulative();
            break;
    }
}

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, logger, particleSizeDistribution_, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_NUMBER_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, validateCumulativeDistribution(), validateProbabilityDensityDistribution(), and WARN.

Referenced by convertPSD2ToPSD(), ParticleHandler::getPSD(), Calibration::setPSD(), Material::setPSD(), setPSDFromCSV(), PSDSelfTest::setupInitialConditions(), and testPSD().

validateCumulativeDistribution()

void PSD::validateCumulativeDistribution

(

)

Validates if a CDF starts with zero and adds up to unity.

Validates if a distribution is cumulative by first checking if the psd vector is empty and that the scaling of probabilities is in the range [0,1]. Also it checks if each consecutive value is higher than the latter value (i.e. assuming piecewise linear CDF: p_i > p_i-1). Further it replaces probabilities if the CDF does not start with zero or does not end with unity (i.e. p_0=0 and p_end=1)

{
    // check whether the distribution is cumulative
    for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
    {
        logger.assert_always(it->probability >= (it - 1)->probability,
                             "psd is not cumulative: radius % %, probabilities % %", (it - 1)->internalVariable,
                             it->internalVariable, (it - 1)->probability, it->probability);
    }
    // set negative probabilities to 0
    for (auto& rp : particleSizeDistribution_)
    {
        if (rp.probability < 0)
        {
            logger(INFO, "negative probability % has been set to 0", rp.probability);
            rp.probability = 0;
        }
    }
    // ensure interval of probabilities to be [0,1]
    for (auto p = 0; p < particleSizeDistribution_.size(); ++p)
        particleSizeDistribution_[p].probability /= particleSizeDistribution_.back().probability;
    // cdf needs to start with a probability of zero
    if (particleSizeDistribution_[0].probability != 0)
    {
        logger(INFO, "adding a zero at the beginning of the psd");
        particleSizeDistribution_.insert(particleSizeDistribution_.begin(),
                                         {std::max(1e-3 * particleSizeDistribution_[0].internalVariable,
                                                   2 * particleSizeDistribution_[0].internalVariable -
                                                   particleSizeDistribution_[1].internalVariable), 0});
    }
    // cdf needs to end with a probability of one
    if (particleSizeDistribution_.back().probability < 1)
    {
        logger(INFO, "adding a one at the end of the psd");
        particleSizeDistribution_.push_back({particleSizeDistribution_.back().internalVariable, 1});
    }
    // cut off equal subsequent values on the end of the cdf to remove zero size classes.
    for (auto i = particleSizeDistribution_.end() - 1; i != particleSizeDistribution_.begin(); i--)
    {
        if (i->probability == (i - 1)->probability)
        {
            particleSizeDistribution_.erase(particleSizeDistribution_.end() + std::distance(particleSizeDistribution_
                                                                                                    .end(), i));
        }
        else
        {
            break;
        }
    }
}

References i, INFO, logger, p, and particleSizeDistribution_.

Referenced by convertProbabilityDensityToCumulative(), setDistributionLogNormal(), setDistributionNormal(), setDistributionPhiNormal(), setDistributionUniform(), and setPSDFromVector().

validateProbabilityDensityDistribution()

void PSD::validateProbabilityDensityDistribution

(

)

Validates if the integral of the PDF equals to unity.

Validates if a PDF is valid by first checking if the PDF starts with zero (i.e. p_0=0). Further it checks if the integral is equal to unity (i.e. assuming piecewise constant PDF: sum(p_i)=1).

{
    // set negative probabilities to 0
    for (auto& rp : particleSizeDistribution_)
    {
        if (rp.probability < 0)
        {
            logger(INFO, "negative probability % has been set to 0", rp.probability);
            rp.probability = 0;
        }
    }
    Mdouble sum = 0;
    // check if the psd starts with zero probability (i.e. p_0=0)
    if (particleSizeDistribution_[0].probability != 0)
    {
        logger(INFO, "adding a zero at the beginning of PDF");
        particleSizeDistribution_.insert(particleSizeDistribution_.begin(),
                                         {std::max(1e-3 * particleSizeDistribution_[0].internalVariable,
                                                   2 * particleSizeDistribution_[0].internalVariable -
                                                   particleSizeDistribution_[1].internalVariable), 0});
    }
    // sum up probabilities to check if it equals to unity (sum(p_i)=1)
    for (auto& it : particleSizeDistribution_)
    {
        sum += it.probability;
    }
    // ensure interval of probabilities to be [0,1] by normalization.
    for (auto& it : particleSizeDistribution_)
        it.probability /= sum;
    // if the sum of probabilities is not equal to unity, sum up the normalized probabilities again
    if (sum - 1 > 1e-6)
    {
        sum = 0;
        for (auto& it: particleSizeDistribution_)
        {
            sum += it.probability;
        }
    }
    logger.assert_debug(sum - 1 < 1e-6, "PDF is not valid: Integral of PDF is not equal to unity");
}

References e(), INFO, logger, and particleSizeDistribution_.

Referenced by convertCumulativeToProbabilityDensity(), and setPSDFromVector().

Friends And Related Function Documentation

operator< [1/2]

bool operator< ( const DistributionElements &  l, const DistributionElements &  r  )

friend

determines if a certain value of the PSD vector is lower than another one. Used for std::lower_bound()

Required to use std::lower_bound for finding when the probability is higher than a certain value.

Returns

TRUE if probability from a vector of type DistributionElements is higher than a certain value from a vector of type DistributionElements and FALSE in the opposite case.

{
    return l.probability < r.probability;
}

operator< [2/2]

bool operator< ( const DistributionElements &  l, Mdouble  r  )

friend

determines if a certain value of the PSD vector is lower than a double.

required to use std::lower_bound for finding when the probability provided as a double is higher than a certain value.

Returns

TRUE if probability as double is higher than a certain value from a DistributionElements vector and FALSE in the opposite case.

{
    return l.probability < prob;
}

operator<<

std::ostream& operator<< ( std::ostream &  os, DistributionElements &  p  )

friend

Writes to output stream.

Writes to output stream. This function is used for restart files.

Returns

a reference to an output stream.

{
    os << p.internalVariable << ' ' << p.probability << ' ';
    return os;
}

operator==

Mdouble operator== ( DistributionElements  l, Mdouble  r  )

friend

Determines if a certain value of the PSD vector is equal to a double.

Required to use std::distance to find the index of the PSD size class in which a particle has to be inserted

Returns

A double which determines the size class (radius) a particle will be inserted to.

{
    return l.internalVariable == r;
}

operator>>

std::istream& operator>> ( std::istream &  is, DistributionElements &  p  )

friend

Reads from input stream.

reads from input stream. This function is used for restart files.

Returns

a reference to an input stream.

{
    is >> p.internalVariable;
    is >> p.probability;
    return is;
}

Member Data Documentation

index_

int PSD::index_

private

Integer which determines the class in which a particle has to be inserted for the manual insertion routine.

Referenced by decrementNParticlesPerClass(), decrementVolumePerClass(), insertManuallyByVolume(), and PSD().

nParticlesPerClass_

std::vector<int> PSD::nParticlesPerClass_

private

Vector of integers which represents the number of inserted particles in each size class. The classes in this vector are defined to contain the particles between size r_i and r_i-1. (e.g. size class 12 consists of particles between size class 12 and 11 of the PDF)

Referenced by decrementNParticlesPerClass(), getInsertedParticleNumber(), and insertManuallyByVolume().

particleSizeDistribution_

std::vector<DistributionElements> PSD::particleSizeDistribution_

private

Vector of the DistributionElements class which stores radii as internalVariable and probabilities of a PSD.

Referenced by convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensity(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), cutoffAndSqueezeCumulativeDistribution(), cutoffCumulativeDistribution(), drawSample(), getMaxRadius(), getMinRadius(), getParticleSizeDistribution(), getQuantileByRadius(), getRadius(), getRadiusByQuantile(), getVolumetricMeanRadius(), insertManuallyByVolume(), printPSD(), PSD(), scaleParticleSize(), scaleParticleSizeAuto(), setDistributionLogNormal(), setDistributionNormal(), setDistributionPhiNormal(), setDistributionUniform(), setParticleSizeDistribution(), setPSDFromVector(), validateCumulativeDistribution(), and validateProbabilityDensityDistribution().

random_

RNG PSD::random_

private

Mercury random number generator object used to draw random numbers from a random initial seed

Referenced by drawSample(), PSD(), and setFixedSeed().

volumePerClass_

std::vector<Mdouble> PSD::volumePerClass_

private

Vector of doubles which stores the volume of inserted particles for each size class. This vector is used in the insertManuallyByVolume() function to check if the volumeAllowed per class is exceeded and thus no further particles should be added to a certain class. The classes in this vector are defined to contain the volume of particles between size r_i and r_i-1. (e.g. size class 12 consists of the particles' volume between size class 12 and 11 of the PDF)

Referenced by decrementVolumePerClass(), and insertManuallyByVolume().

---

The documentation for this class was generated from the following files:

• PSD.h
• PSD.cc