MPS Class Reference

MPS simulation class. More...

#include <mps.hpp>

Collaboration diagram for MPS:

Public Member Functions
	MPS ()=default
	MPS (const Input &input, std::unique_ptr< PressureCalculator::Interface > &&pressureCalculator, std::unique_ptr< SurfaceDetector::Interface > &&surfaceDetector)
void	stepForward ()

Public Attributes
Settings	settings
	Settings for the simulation.
RefValues	refValuesForNumberDensity
	Reference values for the simulation ( \(n^0\), \(\lambda^0\))
RefValues	refValuesForLaplacian
	Reference values for the simulation ( \(n^0\), \(\lambda^0\))
RefValues	refValuesForGradient
	Reference values for the simulation ( \(n^0\), \(\lambda^0\))
Particles	particles
	Particles in the simulation.
Domain	domain {}
	Domain of the simulation.
std::unique_ptr< PressureCalculator::Interface >	pressureCalculator
	Interface for pressure calculation.
double	courant {}
	Maximum courant number among all particles.

Private Member Functions
void	calGravity ()
	calculate gravity term
void	calViscosity (const double &re)
	calculate viscosity term of Navier-Stokes equation
void	moveParticle ()
	move particles in prediction step
void	collision ()
	calculate collision between particles when they are too close
void	calNumberDensity (const double &re)
	calculate number density of each particle
void	setBoundaryCondition ()
	set boundary condition of pressure Poisson equation
bool	isFreeSurface (const Particle &pi)
bool	isParticleDistributionBiased (const Particle &pi)
	check if particle distribution is biased.
void	setMinimumPressure (const double &re)
	set minimum pressure for pressure gradient calculation
void	calPressureGradient (const double &re)
	calculate pressure gradient term
void	moveParticleUsingPressureGradient ()
	move particles in correction step
void	calCourant ()
	calculate Courant number

Private Attributes
NeighborSearcher	neighborSearcher
	Neighbor searcher for neighbor search.
std::unique_ptr< SurfaceDetector::Interface >	surfaceDetector
	Interface for free surface detection.

Detailed Description

MPS simulation class.

Executes the MPS simulation. This class does not handle the simulation process itself, but only the calculation of the MPS method.

Definition at line 25 of file mps.hpp.

Constructor & Destructor Documentation

MPS() [1/2]

MPS::MPS ( )

default

MPS() [2/2]

MPS::MPS	(	const Input &	input,
		std::unique_ptr< PressureCalculator::Interface > &&	pressureCalculator,
		std::unique_ptr< SurfaceDetector::Interface > &&	surfaceDetector )

Definition at line 15 of file mps.cpp.

  {
    this->settings           = input.settings;
    this->domain             = input.settings.domain;
    this->particles          = input.particles;
    this->pressureCalculator = std::move(pressureCalculator);
    this->surfaceDetector    = std::move(surfaceDetector);
    this->neighborSearcher   = NeighborSearcher(input.settings.reMax, input.settings.domain, input.particles.size());
 
    refValuesForNumberDensity = RefValues(settings.dim, settings.particleDistance, settings.re_forNumberDensity);
    refValuesForGradient      = RefValues(settings.dim, settings.particleDistance, settings.re_forGradient);
    refValuesForLaplacian     = RefValues(settings.dim, settings.particleDistance, settings.re_forLaplacian);
}

Here is the call graph for this function:

Member Function Documentation

stepForward()

void MPS::stepForward

(

)

Definition at line 32 of file mps.cpp.

                      {
    neighborSearcher.setNeighbors(particles);
    calGravity();
    calViscosity(settings.re_forLaplacian);
    moveParticle();
 
    neighborSearcher.setNeighbors(particles);
    collision();
 
    neighborSearcher.setNeighbors(particles);
    calNumberDensity(settings.re_forNumberDensity);
    auto pressures = pressureCalculator->calc(particles);
    for (auto& particle : particles) {
        particle.pressure = pressures[particle.id];
    }
 
    setMinimumPressure(settings.re_forGradient);
    calPressureGradient(settings.re_forGradient);
    moveParticleUsingPressureGradient();
 
    // Update pressure again when using EMPS
    if (auto explicitPressureCalculator = dynamic_cast<PressureCalculator::Explicit*>(pressureCalculator.get())) {
        auto pressures = explicitPressureCalculator->calc(particles);
        for (auto& particle : particles) {
            particle.pressure = pressures[particle.id];
        }
    }
 
    calCourant();
}

Here is the call graph for this function:

Here is the caller graph for this function:

calGravity()

void MPS::calGravity ( )

private

calculate gravity term

Definition at line 63 of file mps.cpp.

                     {
#pragma omp parallel for
    for (auto& p : particles) {
        if (p.type == ParticleType::Fluid) {
            p.acceleration += settings.gravity;
 
        } else {
            p.acceleration.setZero();
        }
    }
}

Here is the caller graph for this function:

calViscosity()

void MPS::calViscosity ( const double & re )

private

calculate viscosity term of Navier-Stokes equation

Parameters

re	effective radius \(r_e\)

The viscosity term of the Navier-Stokes equation is calculated as follows:

\[ \nu\langle \nabla^2\mathbf{u}\rangle_i = \nu\frac{2 d}{n^0\lambda^0}\sum_{j\neq i} (\mathbf{u}_j - \mathbf{u}_i) w_{ij} \]

Definition at line 75 of file mps.cpp.

                                       {
    double n0     = refValuesForLaplacian.n0;
    double lambda = refValuesForLaplacian.lambda;
    double a      = (settings.kinematicViscosity) * (2.0 * settings.dim) / (n0 * lambda);
 
#pragma omp parallel for
    for (auto& pi : particles) {
        if (pi.type != ParticleType::Fluid)
            continue;
 
        Eigen::Vector3d viscosityTerm = Eigen::Vector3d::Zero();
 
        for (auto& neighbor : pi.neighbors) {
            auto& pj = particles[neighbor.id];
 
            if (neighbor.distance < settings.re_forLaplacian) {
                double w = weight(neighbor.distance, re);
                viscosityTerm += (pj.velocity - pi.velocity) * w;
            }
        }
 
        viscosityTerm *= a;
        pi.acceleration += viscosityTerm;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

moveParticle()

void MPS::moveParticle ( )

private

move particles in prediction step

The position and velocity of each particle are updated as

\[ \mathbf{u}_i^* = \mathbf{u}_i^k + (\nu\langle \nabla^2\mathbf{u}\rangle^k_i+\mathbf{g})\Delta t. \]

Definition at line 101 of file mps.cpp.

                       {
#pragma omp parallel for
    for (auto& p : particles) {
        if (p.type == ParticleType::Fluid) {
            p.velocity += p.acceleration * settings.dt;
            p.position += p.velocity * settings.dt;
        }
        p.acceleration.setZero();
    }
}

Here is the caller graph for this function:

collision()

void MPS::collision ( )

private

calculate collision between particles when they are too close

Definition at line 112 of file mps.cpp.

                    {
    for (auto& pi : particles) {
        if (pi.type != ParticleType::Fluid)
            continue;
 
        for (auto& neighbor : pi.neighbors) {
            Particle& pj = particles[neighbor.id];
            if (pj.type == ParticleType::Fluid && pj.id >= pi.id)
                continue;
 
            if (neighbor.distance < settings.collisionDistance) {
 
                double invMi = pi.inverseDensity();
                double invMj = pj.inverseDensity();
                double mass  = 1.0 / (invMi + invMj);
 
                Eigen::Vector3d normal = (pj.position - pi.position).normalized();
                double relVel          = (pj.velocity - pi.velocity).dot(normal);
                double impulse         = 0.0;
                if (relVel < 0.0)
                    impulse = -(1 + settings.coefficientOfRestitution) * relVel * mass;
                pi.velocity -= impulse * invMi * normal;
                pj.velocity += impulse * invMj * normal;
 
                double depth           = settings.collisionDistance - neighbor.distance;
                double positionImpulse = depth * mass;
                pi.position -= positionImpulse * invMi * normal;
                pj.position += positionImpulse * invMj * normal;
 
                // cerr << "WARNING: Collision between particles " << pi.id << " and " << pj.id << " occurred."
                // << endl;
            }
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

calNumberDensity()

void MPS::calNumberDensity ( const double & re )

private

calculate number density of each particle

Parameters

re	effective radius \(r_e\)

Definition at line 148 of file mps.cpp.

                                           {
#pragma omp parallel for
    for (auto& pi : particles) {
        pi.numberDensity = 0.0;
 
        if (pi.type == ParticleType::Ghost)
            continue;
 
        for (auto& neighbor : pi.neighbors)
            pi.numberDensity += weight(neighbor.distance, re);
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

setBoundaryCondition()

void MPS::setBoundaryCondition ( )

private

set boundary condition of pressure Poisson equation

Definition at line 161 of file mps.cpp.

                               {
#pragma omp parallel for
    for (auto& pi : particles) {
        if (pi.type == ParticleType::Ghost || pi.type == ParticleType::DummyWall) {
            pi.boundaryCondition = FluidState::Ignored;
 
        } else { // Fluid particles
            if (surfaceDetector->isFreeSurface(particles, pi)) {
                pi.boundaryCondition = FluidState::FreeSurface;
            } else {
                pi.boundaryCondition = FluidState::Inner;
            }
        }
    }
}

isFreeSurface()

bool MPS::isFreeSurface ( const Particle & pi )

private

isParticleDistributionBiased()

bool MPS::isParticleDistributionBiased ( const Particle & pi )

private

check if particle distribution is biased.

Particle number density can be low even in the fluid region, making the free surface detection based on particle number density vulnerable. Even when a specific particle's the particle number density is lower than criteria, it can be considered as inner particle if the particle distribution around the particle is not biased. In this way we can avoid the free surface detection error. This is based on Khayyer et al. (https://doi.org/10.1016/j.apor.2009.06.003)

Parameters

pi	Particle to check

Returns

true particle distribution is biased = free surface

false particle distribution is not biased = not free surface

setMinimumPressure()

void MPS::setMinimumPressure ( const double & re )

private

set minimum pressure for pressure gradient calculation

Parameters

re	effective radius \(r_e\)

Definition at line 177 of file mps.cpp.

                                             {
#pragma omp parallel for
    for (auto& p : particles) {
        p.minimumPressure = p.pressure;
    }
 
    for (auto& pi : particles) {
        if (pi.type == ParticleType::Ghost || pi.type == ParticleType::DummyWall)
            continue;
 
        for (auto& neighbor : pi.neighbors) {
            Particle& pj = particles[neighbor.id];
            if (pj.type == ParticleType::Ghost || pj.type == ParticleType::DummyWall)
                continue;
            if (pj.id > pi.id)
                continue;
 
            if (neighbor.distance < re) {
                pi.minimumPressure = std::min(pi.minimumPressure, pj.pressure);
                pj.minimumPressure = std::min(pj.minimumPressure, pi.pressure);
            }
        }
    }
}

Here is the caller graph for this function:

calPressureGradient()

void MPS::calPressureGradient ( const double & re )

private

calculate pressure gradient term

Parameters

re

The pressure gradient term of the Navier-Stokes equation is calculated as

\[ -\frac{1}{\rho^0}\langle\nabla P\rangle_i = -\frac{1}{\rho^0}\frac{d}{n^0}\sum_{j\neq i} \frac{P_j-P'_i}{\|\mathbf{r}_{ij}\|^2}\mathbf{r}_{ij} w_{ij} \]

where \(P'_i\) is the minimum pressure of the particle \(i\).

Definition at line 202 of file mps.cpp.

                                              {
    double a = settings.dim / refValuesForGradient.n0;
 
#pragma omp parallel for
    for (auto& pi : particles) {
        if (pi.type != ParticleType::Fluid)
            continue;
 
        Eigen::Vector3d grad = Eigen::Vector3d::Zero();
        for (auto& neighbor : pi.neighbors) {
            Particle& pj = particles[neighbor.id];
            if (pj.type == ParticleType::Ghost || pj.type == ParticleType::DummyWall)
                continue;
 
            if (neighbor.distance < re) {
                double w = weight(neighbor.distance, re);
                // double dist2 = pow(neighbor.distance, 2);
                double dist2 = (pj.position - pi.position).squaredNorm();
                double pij   = (pj.pressure - pi.minimumPressure) / dist2;
                grad += (pj.position - pi.position) * pij * w;
            }
        }
        grad *= a;
        pi.acceleration -= grad / pi.density;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

moveParticleUsingPressureGradient()

void MPS::moveParticleUsingPressureGradient ( )

private

move particles in correction step

The position and velocity of each particle are updated as

\[ \mathbf{u}_i^{k+1} = \mathbf{u}_i^* - \frac{1}{\rho^0} \langle\nabla P^{k+1} \rangle_i \Delta t. \]

Definition at line 229 of file mps.cpp.

                                            {
#pragma omp parallel for
    for (auto&& p : particles) {
        if (p.type == ParticleType::Fluid) {
            p.velocity += p.acceleration * settings.dt;
            p.position += p.acceleration * settings.dt * settings.dt;
        }
 
        p.acceleration.setZero();
    }
}

Here is the caller graph for this function:

calCourant()

void MPS::calCourant ( )

private

calculate Courant number

Definition at line 241 of file mps.cpp.

                     {
    courant = 0.0;
 
    for (auto& pi : particles) {
        if (pi.type != ParticleType::Fluid)
            continue;
 
        double iCourant = (pi.velocity.norm() * settings.dt) / settings.particleDistance;
        courant         = std::max(courant, iCourant);
    }
 
    if (courant > settings.cflCondition) {
        cerr << "ERROR: Courant number is larger than CFL condition. Courant = " << courant << endl;
    }
}

Here is the caller graph for this function:

Member Data Documentation

settings

Settings MPS::settings

Settings for the simulation.

Definition at line 27 of file mps.hpp.

refValuesForNumberDensity

RefValues MPS::refValuesForNumberDensity

Reference values for the simulation ( \(n^0\), \(\lambda^0\))

Definition at line 28 of file mps.hpp.

refValuesForLaplacian

RefValues MPS::refValuesForLaplacian

Reference values for the simulation ( \(n^0\), \(\lambda^0\))

Definition at line 29 of file mps.hpp.

refValuesForGradient

RefValues MPS::refValuesForGradient

Reference values for the simulation ( \(n^0\), \(\lambda^0\))

Definition at line 30 of file mps.hpp.

particles

Particles MPS::particles

Particles in the simulation.

Definition at line 31 of file mps.hpp.

domain

Domain MPS::domain {}

Domain of the simulation.

Definition at line 32 of file mps.hpp.

{};                     

pressureCalculator

std::unique_ptr<PressureCalculator::Interface> MPS::pressureCalculator

Interface for pressure calculation.

Definition at line 34 of file mps.hpp.

courant

double MPS::courant {}

Maximum courant number among all particles.

Definition at line 36 of file mps.hpp.

{}; 

neighborSearcher

NeighborSearcher MPS::neighborSearcher

private

Neighbor searcher for neighbor search.

Definition at line 47 of file mps.hpp.

surfaceDetector

std::unique_ptr<SurfaceDetector::Interface> MPS::surfaceDetector

private

Interface for free surface detection.

Definition at line 48 of file mps.hpp.

---

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

• src/mps.hpp
• src/mps.cpp