Class for setting up and solving pressure Poisson equation. More...

#include <pressure_poisson_equation.hpp>

Public Member Functions
	PressurePoissonEquation ()=default

	PressurePoissonEquation (int dimension, double dt, double relaxationCoefficient, double compressibility, double n0_forNumberDensity, double n0_forLaplacian, double lambda0, double reForLaplacian, double reForNumberDensity)

void	setup (const Particles &particles, const DirichletBoundaryCondition &dirichletBoundaryCondition)
	Setup pressure Poisson equation.

std::vector< double >	solve ()
	Solve pressure Poisson equation.

Private Member Functions
void	resetEquation ()

void	setSourceTerm (const Particles &particles, const DirichletBoundaryCondition &dirichletBoundaryCondition)
	Set the source term for the pressure Poisson equation.

void	setMatrixTriplets (const Particles &particles, const DirichletBoundaryCondition &dirichletBoundaryCondition)
	Set the matrix triplets for the pressure Poisson equation.

Private Attributes
int	dimension

double	dt

double	relaxationCoefficient

double	compressibility

double	n0_forNumberDensity

double	n0_forLaplacian

double	lambda0

double	reForLaplacian

double	reForNumberDensity

size_t	particlesCount

std::vector< Eigen::Triplet< double > >	matrixTriplets
	Triplets for coefficient matrix.

Eigen::SparseMatrix< double, Eigen::RowMajor >	coefficientMatrix
	Coefficient matrix for pressure Poisson equation.

Eigen::VectorXd	sourceTerm
	Source term for pressure Poisson equation.

Detailed Description

Class for setting up and solving pressure Poisson equation.

Definition at line 14 of file pressure_poisson_equation.hpp.

Constructor & Destructor Documentation

◆ PressurePoissonEquation() [1/2]

PressureCalculator::PressurePoissonEquation::PressurePoissonEquation ( )

default

◆ PressurePoissonEquation() [2/2]

PressurePoissonEquation::PressurePoissonEquation	(	int	dimension,
		double	dt,
		double	relaxationCoefficient,
		double	compressibility,
		double	n0_forNumberDensity,
		double	n0_forLaplacian,
		double	lambda0,
		double	reForLaplacian,
		double	reForNumberDensity )

Definition at line 9 of file pressure_poisson_equation.cpp.

  {
    this->dimension             = dimension;
    this->dt                    = dt;
    this->relaxationCoefficient = relaxationCoefficient;
    this->compressibility       = compressibility;
    this->n0_forNumberDensity   = n0_forNumberDensity;
    this->n0_forLaplacian       = n0_forLaplacian;
    this->lambda0               = lambda0;
    this->reForLaplacian        = reForLaplacian;
    this->reForNumberDensity    = reForNumberDensity;
}

Member Function Documentation

◆ setup()

void PressurePoissonEquation::setup	(	const Particles &	particles,
		const DirichletBoundaryCondition &	dirichletBoundaryCondition )

Setup pressure Poisson equation.

Parameters

particles	Particles
isPressureUpdateTarget	Function that gets a particle and returns true if the particle is a target for pressure update

Definition at line 31 of file pressure_poisson_equation.cpp.

  {
    this->particlesCount = particles.size();
 
    resetEquation();
    setSourceTerm(particles, dirichletBoundaryCondition);
    setMatrixTriplets(particles, dirichletBoundaryCondition);
    coefficientMatrix.setFromTriplets(matrixTriplets.begin(), matrixTriplets.end());
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ solve()

std::vector< double > PressurePoissonEquation::solve ( )

Solve pressure Poisson equation.

Returns: Calculated pressure. The size of the vector is the same as the number of particles.

Definition at line 42 of file pressure_poisson_equation.cpp.

                                                 {
    using std::cerr;
    using std::endl;
 
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double, Eigen::RowMajor>> solver;
    solver.compute(coefficientMatrix);
    Eigen::VectorXd pressure = solver.solve(sourceTerm);
    if (solver.info() != Eigen::Success) {
        cerr << "Pressure calculation failed." << endl;
        std::exit(-1);
    }
 
    // This conversion is done by giving std::vector the pointers to the first and the last elements of the
    // Eigen::VectorXd. If this type of conversion appears frequently, consider defining a function to convert the
    // vector type.
    std::vector<double> pressureStdVec(pressure.data(), pressure.data() + pressure.size());
    return pressureStdVec;
}

Here is the caller graph for this function:

◆ resetEquation()

void PressurePoissonEquation::resetEquation ( )

private

Definition at line 61 of file pressure_poisson_equation.cpp.

                                            {
    coefficientMatrix.resize(particlesCount, particlesCount);
    sourceTerm.resize(particlesCount);
    matrixTriplets.clear();
}

Here is the caller graph for this function:

◆ setSourceTerm()

void PressurePoissonEquation::setSourceTerm	(	const Particles &	particles,
		const DirichletBoundaryCondition &	dirichletBoundaryCondition )

private

Set the source term for the pressure Poisson equation.

Parameters

particles	Particles
isPressureUpdateTarget	Function that gets a particle and returns true if the particle is a target for pressure update

Definition at line 73 of file pressure_poisson_equation.cpp.

  {
    double n0    = this->n0_forNumberDensity;
    double gamma = this->relaxationCoefficient;
 
#pragma omp parallel for
    for (auto& pi : particles) {
        if (dirichletBoundaryCondition.contains(pi.id)) {
            // When dirichlet boundary condition is set, the source term is the boundary condition value.
            sourceTerm[pi.id] = dirichletBoundaryCondition.value(pi.id);
        } else {
            sourceTerm[pi.id] = gamma * (1.0 / (dt * dt)) * ((pi.numberDensity - n0) / n0);
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ setMatrixTriplets()

void PressurePoissonEquation::setMatrixTriplets	(	const Particles &	particles,
		const DirichletBoundaryCondition &	dirichletBoundaryCondition )

private

Set the matrix triplets for the pressure Poisson equation.

Coefficient matrix is a sparse matrix and made by triplets. Triplets are the elements of the matrix that are not zero and represented by the row index, column index, and the value of the element. This function sets the triplets for the pressure Poisson equation.

Parameters

particles	Particles
isPressureUpdateTarget	Function that gets a particle and returns true if the particle is a target for pressure update

Definition at line 99 of file pressure_poisson_equation.cpp.

  {
    auto a  = 2.0 * dimension / (n0_forLaplacian * lambda0);
    auto re = reForLaplacian;
 
    for (auto& pi : particles) {
        if (dirichletBoundaryCondition.contains(pi.id)) {
            // When Dirichlet boundary conditions are set, only the diagonal term is set to 1 so that the pressure is at
            // the specified value.
            matrixTriplets.emplace_back(pi.id, pi.id, 1.0);
            continue;
        }
 
        double coefficient_ii = 0.0;
        for (auto& neighbor : pi.neighbors) {
            auto& pj = particles[neighbor.id];
            if (pj.boundaryCondition == FluidState::Ignored) {
                continue;
            }
 
            if (neighbor.distance < re) {
                double coefficient_ij = a * weight(neighbor.distance, re) / pi.density;
                matrixTriplets.emplace_back(pi.id, pj.id, -1.0 * coefficient_ij);
                coefficient_ii += coefficient_ij;
            }
        }
        coefficient_ii += (compressibility) / (dt * dt);
        matrixTriplets.emplace_back(pi.id, pi.id, coefficient_ii);
    }
}