Program Listing for File TimeStepPBF.cpp

↰ Return to documentation for file (SPlisHSPlasH/PBF/TimeStepPBF.cpp)

#include "TimeStepPBF.h"
#include "SPlisHSPlasH/TimeManager.h"
#include "TimeIntegration.h"
#include "SPlisHSPlasH/SPHKernels.h"
#include "SimulationDataPBF.h"
#include <iostream>
#include "Utilities/Timing.h"
#include "SPlisHSPlasH/Simulation.h"
#include "Utilities/Counting.h"
#include "SPlisHSPlasH/BoundaryModel_Akinci2012.h"
#include "SPlisHSPlasH/BoundaryModel_Koschier2017.h"
#include "SPlisHSPlasH/BoundaryModel_Bender2019.h"

using namespace SPH;
using namespace std;
using namespace GenParam;

std::string TimeStepPBF::METHOD_NAME = "PBF";
int TimeStepPBF::SOLVER_ITERATIONS = -1;
int TimeStepPBF::MIN_ITERATIONS = -1;
int TimeStepPBF::MAX_ITERATIONS = -1;
int TimeStepPBF::MAX_ERROR = -1;
int TimeStepPBF::VELOCITY_UPDATE_METHOD = -1;
int TimeStepPBF::ENUM_PBF_FIRST_ORDER = -1;
int TimeStepPBF::ENUM_PBF_SECOND_ORDER = -1;


TimeStepPBF::TimeStepPBF() :
    TimeStep()
{
    m_simulationData.init();
    m_iterations = 0;
    m_minIterations = 2;
    m_maxIterations = 100;
    m_maxError = static_cast<Real>(0.01);
    m_velocityUpdateMethod = 0;

    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        model->addField({ "lambda", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getLambda(fluidModelIndex, i); } });
        model->addField({ "deltaX", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDeltaX(fluidModelIndex, i)[0]; } });
    }
}

TimeStepPBF::~TimeStepPBF(void)
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        model->removeFieldByName("lambda");
        model->removeFieldByName("deltaX");
    }
}

void TimeStepPBF::initParameters()
{
    TimeStep::initParameters();

    SOLVER_ITERATIONS = createNumericParameter("iterations", "Iterations", &m_iterations);
    setGroup(SOLVER_ITERATIONS, "Simulation|PBF");
    setDescription(SOLVER_ITERATIONS, "Iterations required by the pressure solver.");
    getParameter(SOLVER_ITERATIONS)->setReadOnly(true);

    MIN_ITERATIONS = createNumericParameter("minIterations", "Min. iterations", &m_minIterations);
    setGroup(MIN_ITERATIONS, "Simulation|PBF");
    setDescription(MIN_ITERATIONS, "Minimal number of iterations of the pressure solver.");
    static_cast<NumericParameter<unsigned int>*>(getParameter(MIN_ITERATIONS))->setMinValue(0);

    MAX_ITERATIONS = createNumericParameter("maxIterations", "Max. iterations", &m_maxIterations);
    setGroup(MAX_ITERATIONS, "Simulation|PBF");
    setDescription(MAX_ITERATIONS, "Maximal number of iterations of the pressure solver.");
    static_cast<NumericParameter<unsigned int>*>(getParameter(MAX_ITERATIONS))->setMinValue(1);

    MAX_ERROR = createNumericParameter("maxError", "Max. density error(%)", &m_maxError);
    setGroup(MAX_ERROR, "Simulation|PBF");
    setDescription(MAX_ERROR, "Maximal density error (%).");
    static_cast<RealParameter*>(getParameter(MAX_ERROR))->setMinValue(static_cast<Real>(1e-6));

    VELOCITY_UPDATE_METHOD = createEnumParameter("velocityUpdateMethod", "Velocity update method", &m_velocityUpdateMethod);
    setGroup(VELOCITY_UPDATE_METHOD, "Simulation|PBF");
    setDescription(VELOCITY_UPDATE_METHOD, "Method for the velocity integration.");
    EnumParameter *enumParam = static_cast<EnumParameter*>(getParameter(VELOCITY_UPDATE_METHOD));
    enumParam->addEnumValue("First Order Update", ENUM_PBF_FIRST_ORDER);
    enumParam->addEnumValue("Second Order Update", ENUM_PBF_SECOND_ORDER);
}

void TimeStepPBF::step()
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    TimeManager *tm = TimeManager::getCurrent ();
    const Real h = tm->getTimeStepSize();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        clearAccelerations(fluidModelIndex);

        // Time integration
        #pragma omp parallel default(shared)
        {
            #pragma omp for schedule(static)
            for (int i = 0; i < (int) model->numActiveParticles(); i++)
            {
                m_simulationData.getLastPosition(fluidModelIndex, i) = m_simulationData.getOldPosition(fluidModelIndex, i);
                m_simulationData.getOldPosition(fluidModelIndex, i) = model->getPosition(i);
                if (model->getParticleState(i) == ParticleState::Active)
                    TimeIntegration::semiImplicitEuler(h, model->getMass(i), model->getPosition(i), model->getVelocity(i), model->getAcceleration(i));
            }
        }
    }

    // Perform neighborhood search
    sim->performNeighborhoodSearch();

#ifdef USE_PERFORMANCE_OPTIMIZATION
    precomputeValues();
#endif

    // Solve density constraint
    START_TIMING("pressureSolve");
    pressureSolve();
    STOP_TIMING_AVG;

    // Update velocities
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        if (m_velocityUpdateMethod == ENUM_PBF_FIRST_ORDER)
        {
            #pragma omp parallel default(shared)
            {
                #pragma omp for schedule(static)
                for (int i = 0; i < (int)model->numActiveParticles(); i++)
                {
                    if (model->getParticleState(i) == ParticleState::Active)
                        TimeIntegration::velocityUpdateFirstOrder(h, model->getMass(i), model->getPosition(i), m_simulationData.getOldPosition(fluidModelIndex, i), model->getVelocity(i));
                    model->getAcceleration(i).setZero();
                }
            }
        }
        else
        {
            #pragma omp parallel default(shared)
            {
                #pragma omp for schedule(static)
                for (int i = 0; i < (int)model->numActiveParticles(); i++)
                {
                    if (model->getParticleState(i) == ParticleState::Active)
                        TimeIntegration::velocityUpdateSecondOrder(h, model->getMass(i), model->getPosition(i), m_simulationData.getOldPosition(fluidModelIndex, i), m_simulationData.getLastPosition(fluidModelIndex, i), model->getVelocity(i));
                    model->getAcceleration(i).setZero();
                }
            }
        }
    }

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
        computeDensities(fluidModelIndex);
    sim->computeNonPressureForces();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        #pragma omp parallel default(shared)
        {
            #pragma omp for schedule(static)
            for (int i = 0; i < (int)model->numActiveParticles(); i++)
            {
                if (model->getParticleState(i) == ParticleState::Active)
                    model->getVelocity(i) += h * model->getAcceleration(i);
            }
        }
    }

    sim->updateTimeStepSize();

    sim->emitParticles();
    sim->animateParticles();

    // Compute new time
    tm->setTime (tm->getTime () + h);
}


void TimeStepPBF::reset()
{
    TimeStep::reset();
    m_simulationData.reset();
    m_iterations = 0;
}


void TimeStepPBF::pressureSolve()
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nFluids = sim->numberOfFluidModels();

    Real avg_density_err = 0;
    m_iterations = 0;
    bool chk = false;

    while ((!chk || (m_iterations < m_minIterations)) && (m_iterations < m_maxIterations))
    {
        chk = true;
        for (unsigned int i = 0; i < nFluids; i++)
        {
            FluidModel *model = sim->getFluidModel(i);
            const Real density0 = model->getDensity0();

            avg_density_err = 0.0;
            pressureSolveIteration(i, avg_density_err);

            // Maximal allowed density fluctuation
            const Real eta = m_maxError * static_cast<Real>(0.01) * density0;  // maxError is given in percent
            chk = chk && (avg_density_err <= eta);
        }
        m_iterations++;
    }
    INCREASE_COUNTER("PBF - iterations", static_cast<Real>(m_iterations));
}

void TimeStepPBF::pressureSolveIteration(const unsigned int fluidModelIndex, Real &avg_density_err)
{
    Simulation *sim = Simulation::getCurrent();
    FluidModel *model = sim->getFluidModel(fluidModelIndex);
    const unsigned int numParticles = model->numActiveParticles();
    if (numParticles == 0)
        return;

    const Real invH = static_cast<Real>(1.0) / TimeManager::getCurrent()->getTimeStepSize();
    const Real invH2 = invH*invH;
    const unsigned int nFluids = sim->numberOfFluidModels();
    const unsigned int nBoundaries = sim->numberOfBoundaryModels();
    const Real eps = static_cast<Real>(1.0e-6);

    const Real density0 = model->getDensity0();

    #pragma omp parallel default(shared)
    {
        #pragma omp for schedule(static)
        for (int i = 0; i < (int) numParticles; i++)
        {
            Real &density = model->getDensity(i);
            m_simulationData.getLambda(fluidModelIndex, i) = 0.0;

            // Compute current density for particle i
            density = model->getVolume(i) * sim->W_zero();
            const Vector3r &xi = model->getPosition(i);

            if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
                computeVolumeAndBoundaryX(fluidModelIndex, i, xi);
            else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
                computeDensityAndGradient(fluidModelIndex, i, xi);

            // Fluid
            forall_fluid_neighbors(
                density += fm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
            );

            // Boundary
            if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
            {
                forall_boundary_neighbors(
                    density += bm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
                );
            }
            else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
            {
                forall_density_maps(
                    density += rho;
                );
            }
            else // Bender2019
            {
                forall_volume_maps(
                    density += Vj * sim->W(xi - xj);
                );
            }

            const Real density_err = density0 * (max(density, static_cast<Real>(1.0)) - static_cast<Real>(1.0));
            #pragma omp atomic
            avg_density_err += density_err;

            // Evaluate constraint function
            const Real C = std::max(density - static_cast<Real>(1.0), static_cast<Real>(0.0));          // clamp to prevent particle clumping at surface

            if (C != 0.0)
            {
                // Compute gradients dC/dx_j
                Real sum_grad_C2 = 0.0;
                Vector3r gradC_i(0.0, 0.0, 0.0);

                // Fluid
                forall_fluid_neighbors(
                    const Vector3r gradC_j = -fm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
                    sum_grad_C2 += gradC_j.squaredNorm();
                    gradC_i -= gradC_j;
                );

                // Boundary
                if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
                {
                    forall_boundary_neighbors(
                        // Boundary: Akinci2012
                        const Vector3r gradC_j = -bm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
                        gradC_i -= gradC_j;
                    );
                }
                else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
                {
                    forall_density_maps(
                        gradC_i -= gradRho;
                    );
                }
                else   // Bender2019
                {
                    forall_volume_maps(
                        const Vector3r gradC_j = -Vj * sim->gradW(xi - xj);
                        gradC_i -= gradC_j;
                    );
                }

                sum_grad_C2 += gradC_i.squaredNorm();

                // Compute lambda
                Real &lambda = m_simulationData.getLambda(fluidModelIndex, i);
                lambda = -C / (sum_grad_C2 + eps);
            }
            else
                m_simulationData.getLambda(fluidModelIndex, i) = 0.0;
        }
    }

    avg_density_err /= numParticles;

    #pragma omp parallel default(shared)
    {
        #pragma omp for schedule(static)
        for (int i = 0; i < (int)numParticles; i++)
        {
            Vector3r &corr = m_simulationData.getDeltaX(fluidModelIndex, i);

            // Compute position correction
            corr.setZero();

            const Vector3r &xi = model->getPosition(i);

            // Fluid
            forall_fluid_neighbors(
                const Vector3r gradC_j = -fm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
                corr -= (m_simulationData.getLambda(fluidModelIndex, i) + (fm_neighbor->getDensity0() / density0) * m_simulationData.getLambda(pid, neighborIndex)) * gradC_j;
            );

            // Boundary
            if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
            {
                forall_boundary_neighbors(
                    // Boundary: Akinci2012
                    const Vector3r gradC_j = -bm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
                    const Vector3r dx = m_simulationData.getLambda(fluidModelIndex, i) * gradC_j;
                    corr -= dx;

                    bm_neighbor->addForce(xj, model->getMass(i) * dx * invH2);
                );
            }
            else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
            {
                forall_density_maps(
                    const Vector3r dx = m_simulationData.getLambda(fluidModelIndex, i) * gradRho;
                    corr -= dx;
                    bm_neighbor->addForce(xj, model->getMass(i) * dx * invH2);
                );
            }
            else   // Bender2019
            {
                forall_volume_maps(
                    const Vector3r gradC_j = -Vj * sim->gradW(xi - xj);
                    const Vector3r dx = m_simulationData.getLambda(fluidModelIndex, i) * gradC_j;
                    corr -= dx;
                    bm_neighbor->addForce(xj, model->getMass(i) * dx * invH2);
                );
            }
        }
    }

    #pragma omp parallel default(shared)
    {
        #pragma omp for schedule(static)
        for (int i = 0; i < (int)numParticles; i++)
        {
            if (model->getParticleState(i) == ParticleState::Active)
                model->getPosition(i) += m_simulationData.getDeltaX(fluidModelIndex, i);
        }
    }
}

void TimeStepPBF::performNeighborhoodSearchSort()
{
    m_simulationData.performNeighborhoodSearchSort();
}

void TimeStepPBF::emittedParticles(FluidModel *model, const unsigned int startIndex)
{
    m_simulationData.emittedParticles(model, startIndex);
}

void TimeStepPBF::resize()
{
    m_simulationData.init();
}