PSD Class Reference

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

#include <PSD.h>

Classes
class	RadiusAndProbability
	Class which stores radii and probabilities of a PSD. This class should be used as a vector<PSD::RadiusAndProbability>. More...

Public Types
enum	TYPE {   TYPE::CUMULATIVE_NUMBER_DISTRIBUTION, TYPE::CUMULATIVE_LENGTH_DISTRIBUTION, TYPE::CUMULATIVE_VOLUME_DISTRIBUTION, TYPE::CUMULATIVE_AREA_DISTRIBUTION,   TYPE::PROBABILITYDENSITY_NUMBER_DISTRIBUTION, TYPE::PROBABILITYDENSITY_LENGTH_DISTRIBUTION, TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION, TYPE::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 (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...
MERCURY_DEPRECATED void	setPSDFromVector (std::vector< RadiusAndProbability > psd, TYPE PSDType)
	Deprecated version of reading in PSDs from a vector. 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	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	convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution ()
	convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION. More...
void	convertCumulativeToCumulativeNumberDistribution (TYPE CDFType)
	convert any other CDF to a CUMULATIVE_NUMBER_DISTRIBUTION. More...
void	cutoffCumulativeNumber (Mdouble quantileMin, Mdouble quantileMax, Mdouble minPolydispersity=0.1)
	cutoff the PSD at given quantiles. More...
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...
Mdouble	getMinRadius () const
	Get smallest radius of the PSD. More...
Mdouble	getMaxRadius () const
	Get largest radius of the PSD. More...
std::vector< RadiusAndProbability >	getParticleSizeDistribution () const
	Get the PSD vector. More...
void	setParticleSizeDistribution (std::vector< RadiusAndProbability >)
	set the PSD by a suitable vector. More...
int	getInsertedParticleNumber () const
	Get the number of particles already inserted into the simulation. 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	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	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...
void	computeRawMomenta ()
	compute raw momenta of the user defined PSD. More...
void	computeCentralMomenta ()
	compute central momenta of the user defined PSD. More...
void	computeStandardisedMomenta ()
	compute standardised momenta of the user defined PSD. More...
std::array< Mdouble, 6 >	getMomenta () const
	get momenta of the user defined PSD. More...

Private Attributes
std::vector< RadiusAndProbability >	particleSizeDistribution_
std::array< Mdouble, 6 >	momenta_ {}
std::vector< int >	nParticlesPerClass_
std::vector< Mdouble >	volumePerClass_

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

Definition at line 62 of file PSD.h.

Member Enumeration Documentation

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

Definition at line 71 of file PSD.h.

    {
        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::PSD

(

)

Constructor; sets everything to 0 or default.

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

Definition at line 32 of file PSD.cc.

References momenta_.

Referenced by copy().

{
    for (auto& momenta: momenta_)
        momenta = 0;
}

PSD::PSD

(

const PSD &

other

)

Copy constructor with deep copy.

Copy constructor

Definition at line 41 of file PSD.cc.

References momenta_, and particleSizeDistribution_.

{
    particleSizeDistribution_ = other.particleSizeDistribution_;
    momenta_ = other.momenta_;
}

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

void PSD::computeCentralMomenta

(

)

compute central momenta of the user defined PSD.

Compute the central momenta of the PSD from their respective raw momenta.

Definition at line 876 of file PSD.cc.

References computeRawMomenta(), and momenta_.

Referenced by computeStandardisedMomenta().

{
    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;
}

void PSD::computeRawMomenta

(

)

compute raw momenta of the user defined PSD.

Compute the raw momenta of the inserted PSD by converting it to a PDF, calculating the moments m_i according to \( m_i = \sum_{j=0}^n \) 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 Wikipedia for details.

Definition at line 856 of file PSD.cc.

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityToCumulative(), getInsertedParticleNumber(), and momenta_.

Referenced by computeCentralMomenta().

{
    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.radius, im) * it.probability;
        }
    }
    // zeroth-moment equals particle number for a PROBABILITYDENSITY_NUMBER_DISTRIBUTION
    momenta_[0] *= getInsertedParticleNumber();
    convertProbabilityDensityToCumulative();
}

void PSD::computeStandardisedMomenta

(

)

compute standardised momenta of the user defined PSD.

Compute the standardised momenta of the PSD from their respective central momenta.

Definition at line 890 of file PSD.cc.

References computeCentralMomenta(), and momenta_.

{
    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;
}

void PSD::convertCumulativeToCumulativeNumberDistribution

(

TYPE

CDFType

)

convert any other CDF to a CUMULATIVE_NUMBER_DISTRIBUTION.

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

Parameters

[in]

PDFType

Type of the PDF: PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_AREA_DISTRIBUTION or PROBABILITYDENSITY_VOLUME_DISTRIBUTION, Where L = Length, A = Area and V = Volume.

Todo:

TP: NOT WORKING! If anyone knows how to do it feel free to add

Definition at line 609 of file PSD.cc.

{
    Mdouble sum = 0;
    switch (CDFType)
    {
        default:
            logger(ERROR, "Wrong CDFType");
            break;
        case TYPE::CUMULATIVE_LENGTH_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                // add conversion here
            
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
        case TYPE::CUMULATIVE_AREA_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                // add conversion here
            
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
        case TYPE::CUMULATIVE_VOLUME_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                // add conversion here
            
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
    }
}

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).

Definition at line 502 of file PSD.cc.

References particleSizeDistribution_, and validateProbabilityDensityDistribution().

Referenced by computeRawMomenta(), getVolumeDx(), getVolumetricMeanRadius(), insertManuallyByVolume(), setPSDFromCSV(), and setPSDFromVector().

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

void PSD::convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution

(

)

convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION.

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

Definition at line 584 of file PSD.cc.

References particleSizeDistribution_, and mathsFunc::square().

Referenced by getVolumeDx(), getVolumetricMeanRadius(), and insertManuallyByVolume().

{
    Mdouble sum = 0;
    for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
    {
        it->probability *=
                0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
        // old conversion
        //p.probability /= cubic(p.radius);
        // sum up probabilities
        sum += it->probability;
    }
    // normalize
    for (auto& p : particleSizeDistribution_)
    {
        p.probability /= sum;
    }
}

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).

Definition at line 487 of file PSD.cc.

References particleSizeDistribution_, and validateCumulativeDistribution().

Referenced by computeRawMomenta(), getVolumeDx(), insertManuallyByVolume(), setPSDFromCSV(), and setPSDFromVector().

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

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, Where L = Length, A = Area and V = Volume.

Definition at line 522 of file PSD.cc.

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

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

{
    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->radius + (it - 1)->radius);
                // old conversion
                //p.probability /= p.radius;
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
        case TYPE::PROBABILITYDENSITY_AREA_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability /=
                        1.0 / 3.0 * (square(it->radius) + it->radius * (it - 1)->radius + square((it - 1)->radius));
                // old conversion
                //p.probability /= square(p.radius);
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
        case TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
            {
                it->probability /=
                        0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
                // old conversion
                //p.probability /= cubic(p.radius);
                // sum up probabilities
                sum += it->probability;
            }
            // normalize
            for (auto& p : particleSizeDistribution_)
            {
                p.probability /= sum;
            }
            break;
    }
}

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

Definition at line 59 of file PSD.cc.

References PSD().

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

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.

Definition at line 789 of file PSD.cc.

References cutoffCumulativeNumber(), getSizeRatio(), logger, and WARN.

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

void PSD::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.

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]	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 of polydispersity ([0,1]) at the base of the CDF.

Definition at line 702 of file PSD.cc.

References cutoffCumulativeNumber(), and getNumberDx().

{
    Mdouble r50 = 0.5 * PSD::getNumberDx(50);
    // cut off
    cutoffCumulativeNumber(quantileMin, quantileMax, minPolydispersity);
    // squeeze psd
    for (auto& p: particleSizeDistribution_)
    {
        p.radius = r50 + (p.radius - r50) * squeeze;
    }
}

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

cutoff the PSD at given quantiles.

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

Parameters

[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 of polydispersity ([0,1]) at the base of the CDF.

Definition at line 668 of file PSD.cc.

References getRadiusByQuantile(), getSizeRatio(), INFO, and logger.

Referenced by cutHighSizeRatio(), and cutoffAndSqueezeCumulative().

{
    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().radius <= radiusMinCut)
    {
        particleSizeDistribution_.erase(particleSizeDistribution_.begin());
    }
    particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMinCut, quantileMin});
    particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMin, 0});
    // cut off max
    while (particleSizeDistribution_.back().radius >= radiusMax)
    {
        particleSizeDistribution_.pop_back();
    }
    Mdouble radiusMaxCut = std::max(radiusMax - (1 - minPolydispersity) * (radiusMax - particleSizeDistribution_.back()
            .radius), radiusMin);
    particleSizeDistribution_.push_back({radiusMaxCut, quantileMax});
    particleSizeDistribution_.push_back({radiusMax, 1});
    logger(INFO, "PSD was cut successfully and now has a size ratio of %", getSizeRatio());
}

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.

Todo:

TP: We should add a variable seed. Maybe make it definable by the user?

Definition at line 109 of file PSD.cc.

References particleSizeDistribution_.

{
    // draw a number between 0 and 1, uniformly distributed
    static std::mt19937 gen(0);
    static std::uniform_real_distribution<Mdouble> dist(0, 1);
    Mdouble prob = dist(gen);
    // 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->radius;
    Mdouble rMax = high->radius;
    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;
}

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.

Definition at line 840 of file PSD.cc.

References nParticlesPerClass_.

Referenced by computeRawMomenta().

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

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.

Definition at line 812 of file PSD.cc.

Referenced by getSizeRatio().

{
        return particleSizeDistribution_.back().radius;
}

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.

Definition at line 803 of file PSD.cc.

Referenced by getSizeRatio().

{
    return particleSizeDistribution_[0].radius;
}

std::array< Mdouble, 6 > PSD::getMomenta

(

)

const

get momenta of the user defined PSD.

Get momenta of the user defined particleSizeDistributionVector_ vector.

Returns

Array of momenta corresponding to the first six moments of the user defined PSD.

Definition at line 903 of file PSD.cc.

References momenta_.

{
    return momenta_;
}

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 which determines the obtained diameter as a percentile of the PSD.

Definition at line 720 of file PSD.cc.

References getRadiusByQuantile().

Referenced by cutoffAndSqueezeCumulative().

{
    return 2.0 * getRadiusByQuantile(x / 100);
}

std::vector< PSD::RadiusAndProbability > 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.

Definition at line 821 of file PSD.cc.

References particleSizeDistribution_.

{
    return particleSizeDistribution_;
}

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.

Definition at line 744 of file PSD.cc.

References logger.

Referenced by cutoffCumulativeNumber(), getNumberDx(), and getVolumeDx().

{
    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->radius;
    else
        return low->radius + (high->radius - low->radius) * (quantile - low->probability) /
                             (high->probability - low->probability);
}

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.

Definition at line 778 of file PSD.cc.

References getMaxRadius(), and getMinRadius().

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

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

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 which determines the obtained diameter as a percentile of the PSD.

Definition at line 730 of file PSD.cc.

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), convertProbabilityDensityToCumulative(), and getRadiusByQuantile().

{
    PSD psd = *this;
    psd.convertCumulativeToProbabilityDensity();
    psd.convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution();
    psd.convertProbabilityDensityToCumulative();
    return 2.0 * psd.getRadiusByQuantile(x / 100);
}

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.

Definition at line 762 of file PSD.cc.

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), and mathsFunc::cubic().

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

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.

Definition at line 137 of file PSD.cc.

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), nParticlesPerClass_, particleSizeDistribution_, constants::pi, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, and volumePerClass_.

{
    // 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->radius;
        Mdouble rMax = high->radius;
        Mdouble pMin = low->probability;
        Mdouble pMax = high->probability;
        // interpolate linearly between [low.radius,high.radius] (assumption: CDF is piecewise linear)
        int index = std::distance(particleSizeDistribution_.begin(), high);
        Mdouble a = (prob - pMin) / (pMax - pMin);
        Mdouble rad = a * rMin + (1 - a) * rMax;
        // Compare inserted volume against volumeAllowed
        Mdouble radVolume = 4.0 / 3.0 * constants::pi * rad * rad * rad;
        convertCumulativeToProbabilityDensity();
        convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution();
        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;
}

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.

Definition at line 71 of file PSD.cc.

References INFO, logger, and particleSizeDistribution_.

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

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

Definition at line 314 of file PSD.cc.

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

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

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

Definition at line 290 of file PSD.cc.

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

Referenced by CubeInsertionBoundary::set(), Chute::setupInitialConditions(), and ChuteWithHopper::setupInitialConditions().

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

void PSD::setParticleSizeDistribution

(

std::vector< RadiusAndProbability >

particleSizeDistribution

)

set the PSD by a suitable vector.

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

Parameters

[in]

particleSizeDistribution

vector containing the radii and probabilities specifying the PSD.

Definition at line 830 of file PSD.cc.

{
    particleSizeDistribution_ = particleSizeDistribution;
}

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.

Definition at line 418 of file PSD.cc.

References convertCumulativeToProbabilityDensity(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), CUMULATIVE_AREA_DISTRIBUTION, CUMULATIVE_LENGTH_DISTRIBUTION, CUMULATIVE_VOLUME_DISTRIBUTION, csvReader::getFirstColumn(), csvReader::getSecondColumn(), getSizeRatio(), constants::i, INFO, logger, particleSizeDistribution_, PROBABILITYDENSITY_AREA_DISTRIBUTION, PROBABILITYDENSITY_LENGTH_DISTRIBUTION, PROBABILITYDENSITY_NUMBER_DISTRIBUTION, PROBABILITYDENSITY_VOLUME_DISTRIBUTION, csvReader::read(), csvReader::setHeader(), validateCumulativeDistribution(), and validateProbabilityDensityDistribution().

{
    // 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
    for (int i = 0; i < radii.size(); ++i)
    {
        particleSizeDistribution_.push_back({radii[i], probabilities[i]});
    }
    logger.assert_always(!particleSizeDistribution_.empty(), "PSD cannot be empty");
    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;
    }
    logger(INFO, "A PSD with size ratio % is now ready to be set", getSizeRatio());
}

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

Deprecated version of reading in PSDs from a vector.

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.

Deprecated:

This is the old way of inserting PSDs. In the future use setPSDFromCSV().

Definition at line 352 of file PSD.cc.

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

{
    particleSizeDistribution_ = psdVector;
    switch (PSDType)
    {
        // default case is a cumulative number distribution (CND)
        default:
            validateCumulativeDistribution();
            break;
        case TYPE::CUMULATIVE_VOLUME_DISTRIBUTION:
            validateCumulativeDistribution();
            convertCumulativeToProbabilityDensity();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(
                    TYPE::PROBABILITYDENSITY_VOLUME_DISTRIBUTION);
            convertProbabilityDensityToCumulative();
            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::PROBABILITYDENSITY_VOLUME_DISTRIBUTION:
            validateProbabilityDensityDistribution();
            convertProbabilityDensityToProbabilityDensityNumberDistribution(PSDType);
            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;
    }
}

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)

Definition at line 210 of file PSD.cc.

References constants::i, INFO, logger, and particleSizeDistribution_.

Referenced by convertProbabilityDensityToCumulative(), setDistributionNormal(), setDistributionUniform(), setPSDFromCSV(), and setPSDFromVector().

{
    // ensure interval of probabilities to be [0,1]
    for (auto p = 0; p < particleSizeDistribution_.size(); ++p)
        particleSizeDistribution_[p].probability =
                particleSizeDistribution_[p].probability / particleSizeDistribution_.back().probability;
    // 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", true);
    }
    // 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].radius,
                                                   2 * particleSizeDistribution_[0].radius -
                                                   particleSizeDistribution_[1].radius), 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().radius, 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));
        }
    }
}

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).

Definition at line 251 of file PSD.cc.

References INFO, logger, and particleSizeDistribution_.

Referenced by convertCumulativeToProbabilityDensity(), setPSDFromCSV(), and setPSDFromVector().

{
    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].radius,
                                                   2 * particleSizeDistribution_[0].radius -
                                                   particleSizeDistribution_[1].radius), 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");
}

Friends And Related Function Documentation

bool operator< ( const PSD::RadiusAndProbability &  l, const PSD::RadiusAndProbability &  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 PSD::RadiusAndProbability is higher than a certain value from a vector of type PSD::RadiusAndProbability and FALSE in the opposite case.

Definition at line 913 of file PSD.cc.

{
    return l.probability < r.probability;
}

bool operator< ( const PSD::RadiusAndProbability &  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 PSD::RadiusAndProbability vector and FALSE in the opposite case.

Definition at line 924 of file PSD.cc.

{
    return l.probability < prob;
}

std::ostream& operator<< ( std::ostream &  os, PSD::RadiusAndProbability &  p  )

friend

Writes to output stream.

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

Returns

a reference to an output stream.

Definition at line 953 of file PSD.cc.

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

Mdouble operator== ( PSD::RadiusAndProbability  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.

Definition at line 933 of file PSD.cc.

{
    return l.radius == r;
}

std::istream& operator>> ( std::istream &  is, PSD::RadiusAndProbability &  p  )

friend

Reads from input stream.

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

Returns

a reference to an input stream.

Definition at line 942 of file PSD.cc.

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

Member Data Documentation

std::array<Mdouble, 6> PSD::momenta_ {}

private

Array of doubles which stores the moments of a user defined discrete PROBABILITYDENSITY_NUMBER_DISTRIBUTION.

Todo:

TW can we make this a local variable instead of a class variable?

Definition at line 307 of file PSD.h.

Referenced by computeCentralMomenta(), computeRawMomenta(), computeStandardisedMomenta(), getMomenta(), and PSD().

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)

Definition at line 314 of file PSD.h.

Referenced by getInsertedParticleNumber(), and insertManuallyByVolume().

std::vector<RadiusAndProbability> PSD::particleSizeDistribution_

private

Vector of the PSD::RadiusAndProbability class which stores radii and probabilities of the PSD.

Definition at line 301 of file PSD.h.

Referenced by convertCumulativeToProbabilityDensity(), convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), convertProbabilityDensityToCumulative(), convertProbabilityDensityToProbabilityDensityNumberDistribution(), drawSample(), getParticleSizeDistribution(), insertManuallyByVolume(), printPSD(), PSD(), setDistributionNormal(), setDistributionUniform(), setPSDFromCSV(), setPSDFromVector(), validateCumulativeDistribution(), and validateProbabilityDensityDistribution().

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)

Definition at line 323 of file PSD.h.

Referenced by insertManuallyByVolume().

---

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

• Math/PSD.h
• Math/PSD.cc