.. _program_listing_file_SPlisHSPlasH_IISPH_TimeStepIISPH.cpp:

Program Listing for File TimeStepIISPH.cpp
==========================================

|exhale_lsh| :ref:`Return to documentation for file <file_SPlisHSPlasH_IISPH_TimeStepIISPH.cpp>` (``SPlisHSPlasH/IISPH/TimeStepIISPH.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "TimeStepIISPH.h"
   #include "SPlisHSPlasH/TimeManager.h"
   #include "SPlisHSPlasH/SPHKernels.h"
   #include "SimulationDataIISPH.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;
   
   int TimeStepIISPH::SOLVER_ITERATIONS = -1;
   int TimeStepIISPH::MIN_ITERATIONS = -1;
   int TimeStepIISPH::MAX_ITERATIONS = -1;
   int TimeStepIISPH::MAX_ERROR = -1;
   
   
   std::string TimeStepIISPH::METHOD_NAME = "IISPH";
   
   TimeStepIISPH::TimeStepIISPH() :
       TimeStep()
   {
       m_simulationData.init();
       m_iterations = 0;
       m_minIterations = 2;
       m_maxIterations = 100;
       m_maxError = static_cast<Real>(0.01);
   
       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({ "a_ii", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getAii(fluidModelIndex, i); } });
           model->addField({ "d_ii", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDii(fluidModelIndex, i)[0]; } });
           model->addField({ "d_ij*p_j", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDij_pj(fluidModelIndex, i)[0]; } });
           model->addField({ "pressure", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressure(fluidModelIndex, i); }, true });
           //model->addField({ "last pressure", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getLastPressure(fluidModelIndex, i); } });
           model->addField({ "advected density", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDensityAdv(fluidModelIndex, i); } });
           model->addField({ "pressure acceleration", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressureAccel(fluidModelIndex, i)[0]; } });
       }
   }
   
   TimeStepIISPH::~TimeStepIISPH(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("a_ii");
           model->removeFieldByName("d_ii");
           model->removeFieldByName("d_ij*p_j");
           model->removeFieldByName("pressure");
           //model->removeFieldByName("last pressure");
           model->removeFieldByName("advected density");
           model->removeFieldByName("pressure acceleration");
       }
   }
   
   void TimeStepIISPH::initParameters()
   {
       TimeStep::initParameters();
   
       SOLVER_ITERATIONS = createNumericParameter("iterations", "Iterations", &m_iterations);
       setGroup(SOLVER_ITERATIONS, "Simulation|IISPH");
       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|IISPH");
       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|IISPH");
       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|IISPH");
       setDescription(MAX_ERROR, "Maximal density error (%).");
       static_cast<RealParameter*>(getParameter(MAX_ERROR))->setMinValue(static_cast<Real>(1e-6));
   }
   
   void TimeStepIISPH::step()
   {
       Simulation *sim = Simulation::getCurrent();
       TimeManager *tm = TimeManager::getCurrent ();
       const Real h = tm->getTimeStepSize();
       const unsigned int nModels = sim->numberOfFluidModels();
   
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
           clearAccelerations(fluidModelIndex);
   
       sim->performNeighborhoodSearch();
   
   #ifdef USE_PERFORMANCE_OPTIMIZATION
       precomputeValues();
   #endif
   
       if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
           computeVolumeAndBoundaryX();
       else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
           computeDensityAndGradient();
   
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
           computeDensities(fluidModelIndex);
   
       sim->computeNonPressureForces();
   
       sim->updateTimeStepSize();
   
       // Solve density constraint 
       START_TIMING("predictAdvection");
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
           predictAdvection(fluidModelIndex);
       STOP_TIMING_AVG;
   
       START_TIMING("pressureSolve");
       pressureSolve();
       STOP_TIMING_AVG;
   
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
           integration(fluidModelIndex);
   
       sim->emitParticles();
       sim->animateParticles();
   
       // Compute new time 
       tm->setTime (tm->getTime () + h);
   }
   
   
   void TimeStepIISPH::reset()
   {
       TimeStep::reset();
       m_simulationData.reset();
       m_iterations = 0;
   }
   
   void TimeStepIISPH::predictAdvection(const unsigned int fluidModelIndex)
   {
       Simulation *sim = Simulation::getCurrent();
       FluidModel *model = sim->getFluidModel(fluidModelIndex);
       const unsigned int numParticles = model->numActiveParticles();
       if (numParticles == 0)
           return;
   
       const unsigned int nFluids = sim->numberOfFluidModels();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
       const Real density0 = model->getDensity0();
   
       const Real h = TimeManager::getCurrent()->getTimeStepSize();
   
       // Predict v_adv
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               Vector3r &vel = model->getVelocity(i);
               const Vector3r &accel = model->getAcceleration(i);
               Vector3r &dii = m_simulationData.getDii(fluidModelIndex, i);
               dii.setZero();
   
               if (model->getParticleState(i) == ParticleState::Active)
                   vel += h * accel;
   
               // Compute d_ii
               const Real density = model->getDensity(i) / density0;
               const Real density2 = density*density;
               const Vector3r &xi = model->getPosition(i);
   
               // Fluid
               forall_fluid_neighbors(
                   dii -= fm_neighbor->getVolume(neighborIndex) / density2 * sim->gradW(xi - xj);
               );
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       dii -= bm_neighbor->getVolume(neighborIndex) / density2 * sim->gradW(xi - xj);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       dii += 1.0 / density2 * gradRho;
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       dii -= Vj / density2 * sim->gradW(xi - xj);
                   );
               }
           }
       }
   
       // Compute rho_adv
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               const Real density = model->getDensity(i) / density0;
               Real &densityAdv = m_simulationData.getDensityAdv(fluidModelIndex, i);
               densityAdv = density;
               const Vector3r &xi = model->getPosition(i);
               const Vector3r &vi = model->getVelocity(i);
   
               // Fluid
               forall_fluid_neighbors(
                   const Vector3r &vj = fm_neighbor->getVelocity(neighborIndex);
                   densityAdv += h*fm_neighbor->getVolume(neighborIndex) * (vi - vj).dot(sim->gradW(xi - xj));
               );
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
                       densityAdv += h*bm_neighbor->getVolume(neighborIndex) * (vi - vj).dot(sim->gradW(xi - xj));
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       Vector3r vj;
                       bm_neighbor->getPointVelocity(xi, vj);
                       densityAdv -= h*(vi - vj).dot(gradRho);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       Vector3r vj;
                       bm_neighbor->getPointVelocity(xj, vj);
                       densityAdv += h*Vj * (vi - vj).dot(sim->gradW(xi - xj));
                   );
               }
   
               const Real &pressure = m_simulationData.getPressure(fluidModelIndex, i);
               Real &lastPressure = m_simulationData.getLastPressure(fluidModelIndex, i);
               lastPressure = static_cast<Real>(0.5)*pressure;
   
               // Compute a_ii
               Real &aii = m_simulationData.getAii(fluidModelIndex, i);
               aii = 0.0;
   
               const Vector3r &dii = m_simulationData.getDii(fluidModelIndex, i);
   
               const Real density2 = density * density;
               const Real dpi = model->getVolume(i) / density2;
   
               // Fluid
               forall_fluid_neighbors(
                   // Compute d_ji
                   const Vector3r kernel = sim->gradW(xi - xj);
                   const Vector3r dji = dpi * kernel;          
                   aii += fm_neighbor->getVolume(neighborIndex) * (dii - dji).dot(kernel);
               );
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       const Vector3r kernel = sim->gradW(xi - xj);
                       const Vector3r dji = dpi * kernel;          
                       aii += bm_neighbor->getVolume(neighborIndex) * (dii - dji).dot(kernel);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       const Vector3r dji = -(1.0 / density2) * gradRho;
                       aii -= (dii - dji).dot(gradRho);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       const Vector3r kernel = sim->gradW(xi - xj);
                       const Vector3r dji = dpi * kernel;
                       aii += Vj * (dii - dji).dot(kernel);
                   );
               }
           }
       }
   }
   
   void TimeStepIISPH::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("IISPH - iterations", static_cast<Real>(m_iterations));
   }
   
   void TimeStepIISPH::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 unsigned int nFluids = sim->numberOfFluidModels();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
   
       const Real density0 = model->getDensity0();
       const Real h = TimeManager::getCurrent()->getTimeStepSize();
       const Real h2 = h*h;
       const Real omega = 0.5;
   
       // Compute dij_pj
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               Vector3r &dij_pj = m_simulationData.getDij_pj(fluidModelIndex, i);
               dij_pj.setZero();
   
               const Vector3r &xi = model->getPosition(i);
   
               // Fluid
               forall_fluid_neighbors(
                   const Real densityj = fm_neighbor->getDensity(neighborIndex) / fm_neighbor->getDensity0();
                   const Real densityj2 = densityj*densityj;
   
                   dij_pj -= fm_neighbor->getVolume(neighborIndex) / densityj2 * m_simulationData.getLastPressure(pid, neighborIndex) * sim->gradW(xi - xj);
               );  
           }
       }
   
       // Compute new pressure
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               Real &pi = m_simulationData.getPressure(fluidModelIndex, i);
               pi = 0.0;
   
               const Real &aii = m_simulationData.getAii(fluidModelIndex, i);
               const Real density = model->getDensity(i) / density0;
               const Vector3r &xi = model->getPosition(i);
   
               const Real density2 = density*density;
               const Real dpi = model->getVolume(i) / density2;
               Real sum = 0.0;
   
               // Fluid
               forall_fluid_neighbors(
                   const Vector3r &d_jk_pk = m_simulationData.getDij_pj(pid, neighborIndex);
   
                   // Compute \sum_{k \neq i} djk*pk
                   // Compute d_ji
                   const Vector3r kernel = sim->gradW(xi - xj);
                   const Vector3r dji = dpi * kernel;
                   const Vector3r d_ji_pi = dji * m_simulationData.getLastPressure(fluidModelIndex, i);
   
                   // \sum ( mj * (\sum dij*pj - djj*pj - \sum_{k \neq i} djk*pk) * sim->gradW)
                   sum += fm_neighbor->getVolume(neighborIndex) * (m_simulationData.getDij_pj(fluidModelIndex, i) - m_simulationData.getDii(pid, neighborIndex)*m_simulationData.getLastPressure(pid, neighborIndex) - (d_jk_pk - d_ji_pi)).dot(kernel);
               );
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       sum += bm_neighbor->getVolume(neighborIndex) * m_simulationData.getDij_pj(fluidModelIndex, i).dot(sim->gradW(xi - xj));
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       sum -= m_simulationData.getDij_pj(fluidModelIndex, i).dot(gradRho);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       sum += Vj * m_simulationData.getDij_pj(fluidModelIndex, i).dot(sim->gradW(xi - xj));
                   );
               }
   
               const Real b = static_cast<Real>(1.0) - m_simulationData.getDensityAdv(fluidModelIndex, i);
   
               const Real &lastPi = m_simulationData.getLastPressure(fluidModelIndex, i);
               const Real denom = aii*h2;
               if (fabs(denom) > 1.0e-9)
                   pi = max((static_cast<Real>(1.0) - omega)*lastPi + omega / denom * (b - h2*sum), static_cast<Real>(0.0));
               else
                   pi = 0.0;
   
               if (pi != 0.0)
               {
                   const Real newDensity = density0 * ((aii*pi + sum)*h2 - b) + density0;
   
                   #pragma omp atomic
                   avg_density_err += newDensity - density0;
               }
           }
       }
   
       for (int i = 0; i < (int)numParticles; i++)
       {
           const Real &pi = m_simulationData.getPressure(fluidModelIndex, i);
           Real &lastPi = m_simulationData.getLastPressure(fluidModelIndex, i);
           lastPi = pi;
       }
   
       avg_density_err /= numParticles;
   }
   
   void TimeStepIISPH::integration(const unsigned int fluidModelIndex)
   {
       Simulation *sim = Simulation::getCurrent();
       FluidModel *model = sim->getFluidModel(fluidModelIndex);
       const unsigned int numParticles = model->numActiveParticles();
       if (numParticles == 0)
           return;
   
       // Compute pressure forces
       computePressureAccels(fluidModelIndex);
   
       Real h = TimeManager::getCurrent()->getTimeStepSize();
   
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int) numParticles; i++)
           {
               if (model->getParticleState(i) == ParticleState::Active)
               {
                   Vector3r &pos = model->getPosition(i);
                   Vector3r &vel = model->getVelocity(i);
                   vel += m_simulationData.getPressureAccel(fluidModelIndex, i) * h;
                   pos += vel * h;
               }
           }
       }
   }
   
   void TimeStepIISPH::computePressureAccels(const unsigned int fluidModelIndex)
   {
       Simulation *sim = Simulation::getCurrent();
       FluidModel *model = sim->getFluidModel(fluidModelIndex);
       const unsigned int numParticles = model->numActiveParticles();
       const unsigned int nFluids = sim->numberOfFluidModels();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
       const Real density0 = model->getDensity0();
   
       // Compute pressure forces
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               const Vector3r &xi = model->getPosition(i);
               const Real density_i = model->getDensity(i) / density0;
   
               Vector3r &ai = m_simulationData.getPressureAccel(fluidModelIndex, i);
               ai.setZero();
   
               const Real density2 = density_i*density_i;
               const Real dpi = m_simulationData.getPressure(fluidModelIndex, i) / density2;
   
               // Fluid
               forall_fluid_neighbors(
                   // Pressure 
                   const Real density_j = fm_neighbor->getDensity(neighborIndex) / fm_neighbor->getDensity0();
                   const Real densityj2 = density_j*density_j;
                   const Real dpj = m_simulationData.getPressure(pid, neighborIndex) / densityj2;
                   ai -= fm_neighbor->getVolume(neighborIndex) * (dpi + fm_neighbor->getDensity0() / density0 * dpj) * sim->gradW(xi - xj);
               );
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       const Vector3r a = bm_neighbor->getVolume(neighborIndex) * dpi * sim->gradW(xi - xj);
                       ai -= a;
                       bm_neighbor->addForce(xj, model->getMass(i) * a);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       const Vector3r a = -dpi * gradRho;
                       ai -= a;
                       bm_neighbor->addForce(xj, model->getMass(i) * a);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       const Vector3r a = Vj * dpi* sim->gradW(xi - xj);
                       ai -= a;
                       bm_neighbor->addForce(xj, model->getMass(i) * a);
                   );
               }
           }
       }
   }
   
   void TimeStepIISPH::performNeighborhoodSearchSort()
   {
       m_simulationData.performNeighborhoodSearchSort();
   }
   
   void TimeStepIISPH::emittedParticles(FluidModel *model, const unsigned int startIndex)
   {
       m_simulationData.emittedParticles(model, startIndex);
   }
   
   void TimeStepIISPH::resize()
   {
       m_simulationData.init();
   }