.. _program_listing_file_SPlisHSPlasH_WCSPH_TimeStepWCSPH.cpp:

Program Listing for File TimeStepWCSPH.cpp
==========================================

|exhale_lsh| :ref:`Return to documentation for file <file_SPlisHSPlasH_WCSPH_TimeStepWCSPH.cpp>` (``SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp``)

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

.. code-block:: cpp

   #include "TimeStepWCSPH.h"
   #include "SPlisHSPlasH/TimeManager.h"
   #include "SPlisHSPlasH/SPHKernels.h"
   #include "SimulationDataWCSPH.h"
   #include <iostream>
   #include "Utilities/Timing.h"
   #include "../Simulation.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 TimeStepWCSPH::METHOD_NAME = "WCSPH";
   int TimeStepWCSPH::STIFFNESS = -1;
   int TimeStepWCSPH::EXPONENT = -1;
   
   
   TimeStepWCSPH::TimeStepWCSPH() :
       TimeStep()
   {
       m_simulationData.init();
   
       m_stiffness = 50.0;
       m_exponent = 7.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({ "pressure", METHOD_NAME, FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressure(fluidModelIndex, i); } });
           model->addField({ "pressure acceleration", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressureAccel(fluidModelIndex, i)[0]; } });
       }
   }
   
   TimeStepWCSPH::~TimeStepWCSPH(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("pressure");
           model->removeFieldByName("pressure acceleration");
       }
   }
   
   void TimeStepWCSPH::initParameters()
   {
       TimeStep::initParameters();
   
       STIFFNESS = createNumericParameter("stiffness", "Stiffness", &m_stiffness);
       setGroup(STIFFNESS, "Simulation|WCSPH");
       setDescription(STIFFNESS, "Stiffness coefficient of EOS.");
       static_cast<RealParameter*>(getParameter(STIFFNESS))->setMinValue(static_cast<Real>(1e-6));
   
       EXPONENT = createNumericParameter("exponent", "Exponent (gamma)", &m_exponent);
       setGroup(EXPONENT, "Simulation|WCSPH");
       setDescription(EXPONENT, "Exponent of EOS.");
       static_cast<RealParameter*>(getParameter(EXPONENT))->setMinValue(static_cast<Real>(1e-6));
   
   }
   
   void TimeStepWCSPH::step()
   {
       Simulation *sim = Simulation::getCurrent();
       const unsigned int nModels = sim->numberOfFluidModels();
       TimeManager *tm = TimeManager::getCurrent ();
       const Real h = tm->getTimeStepSize();
   
       sim->performNeighborhoodSearch();
   
   #ifdef USE_PERFORMANCE_OPTIMIZATION
       precomputeValues();
   #endif
   
       if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
           computeVolumeAndBoundaryX();
       else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
           computeDensityAndGradient();
   
       // Compute accelerations: a(t)
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
       {
           clearAccelerations(fluidModelIndex);
           computeDensities(fluidModelIndex);
       }
       sim->computeNonPressureForces();
   
   
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
       {
           FluidModel *model = sim->getFluidModel(fluidModelIndex);
           const Real density0 = model->getDensity0();
           #pragma omp parallel default(shared)
           {
               #pragma omp for schedule(static)  
               for (int i = 0; i < (int)model->numActiveParticles(); i++)
               {
                   Real &density = model->getDensity(i);
                   density = max(density, density0);
                   m_simulationData.getPressure(fluidModelIndex, i) = m_stiffness * (pow(density / density0, m_exponent) - static_cast<Real>(1.0));
               }
           }
   
           computePressureAccels(fluidModelIndex);
       }
   
       sim->updateTimeStepSize();
   
       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)
                   {
                       Vector3r &pos = model->getPosition(i);
                       Vector3r &vel = model->getVelocity(i);
                       Vector3r &accel = model->getAcceleration(i);
                       accel += m_simulationData.getPressureAccel(fluidModelIndex, i);
                       vel += accel * h;
                       pos += vel * h;
                   }
               }
           }
       }
   
       sim->emitParticles();
       sim->animateParticles();
   
       // Compute new time 
       tm->setTime (tm->getTime () + h);
   }
   
   
   void TimeStepWCSPH::reset()
   {
       TimeStep::reset();
       m_simulationData.reset();
   }
   
   void TimeStepWCSPH::computePressureAccels(const unsigned int fluidModelIndex)
   {
       Simulation *sim = Simulation::getCurrent();
       FluidModel *model = sim->getFluidModel(fluidModelIndex);
       const Real density0 = model->getDensity0();
       const unsigned int numParticles = model->numActiveParticles();
       const unsigned int nFluids = sim->numberOfFluidModels();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
   
       // 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);
   
               Vector3r &ai = m_simulationData.getPressureAccel(fluidModelIndex, i);
               ai.setZero();
   
               const Real dpi = m_simulationData.getPressure(fluidModelIndex, i) / (density_i*density_i);
               // Fluid
               forall_fluid_neighbors(
                   // Pressure 
                   const Real density_j = fm_neighbor->getDensity(neighborIndex) * density0 / fm_neighbor->getDensity0();
                   const Real dpj = m_simulationData.getPressure(pid, neighborIndex) / (density_j*density_j);
                   ai -= density0 * fm_neighbor->getVolume(neighborIndex) * (dpi + dpj) * sim->gradW(xi - xj);
               );
   
               // Boundary
               const Real dpj = m_simulationData.getPressure(fluidModelIndex, i) / (density0*density0);
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       const Vector3r a = density0 * bm_neighbor->getVolume(neighborIndex) * (dpi + dpj)* 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 = -density0 * (dpi + dpj)* gradRho;
                       ai -= a;
                       bm_neighbor->addForce(xj, model->getMass(i) * a);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {               
                   forall_volume_maps(
                       const Vector3r a = density0 * Vj * (dpi + dpj)* sim->gradW(xi - xj);
                       ai -= a;
                       bm_neighbor->addForce(xj, model->getMass(i) * a);
                   );
               }
           }
       }
   }
   
   void TimeStepWCSPH::performNeighborhoodSearchSort()
   {
       m_simulationData.performNeighborhoodSearchSort();
   }
   
   void TimeStepWCSPH::emittedParticles(FluidModel *model, const unsigned int startIndex)
   {
       m_simulationData.emittedParticles(model, startIndex);
   }
   
   void TimeStepWCSPH::resize()
   {
       m_simulationData.init();
   }