.. _program_listing_file_SPlisHSPlasH_Viscosity_Viscosity_Bender2017.cpp:

Program Listing for File Viscosity_Bender2017.cpp
=================================================

|exhale_lsh| :ref:`Return to documentation for file <file_SPlisHSPlasH_Viscosity_Viscosity_Bender2017.cpp>` (``SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp``)

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

.. code-block:: cpp

   #include "Viscosity_Bender2017.h"
   #include "SPlisHSPlasH/TimeManager.h"
   #include "Utilities/Counting.h"
   #include "../Simulation.h"
   #include "SPlisHSPlasH/BoundaryModel_Akinci2012.h"
   #include "SPlisHSPlasH/BoundaryModel_Koschier2017.h"
   #include "SPlisHSPlasH/BoundaryModel_Bender2019.h"
   
   using namespace SPH;
   using namespace GenParam;
   
   int Viscosity_Bender2017::ITERATIONS = -1;
   int Viscosity_Bender2017::MAX_ITERATIONS = -1;
   int Viscosity_Bender2017::MAX_ERROR = -1;
   
   
   Viscosity_Bender2017::Viscosity_Bender2017(FluidModel *model) :
       ViscosityBase(model)
   {
       m_targetStrainRate.resize(model->numParticles(), Vector6r::Zero());
       m_viscosityFactor.resize(model->numParticles(), Matrix6r::Zero());
       m_viscosityLambda.resize(model->numParticles(), Vector6r::Zero());
   
       m_iterations = 0;
       m_maxIter = 50;
       m_maxError = 0.01;
   
       model->addField({ "target strain rate", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_targetStrainRate[i][0]; } });
       model->addField({ "viscosity factor", FieldType::Matrix6, [&](const unsigned int i) -> Real* { return &m_viscosityFactor[i](0,0); } });
       model->addField({ "viscosity lambda", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_viscosityLambda[i][0]; } });
   }
   
   Viscosity_Bender2017::~Viscosity_Bender2017(void)
   {
       m_model->removeFieldByName("target strain rate");
       m_model->removeFieldByName("viscosity factor");
       m_model->removeFieldByName("viscosity lambda");
   
       m_targetStrainRate.clear();
       m_viscosityFactor.clear();
       m_viscosityLambda.clear();
   }
   
   void Viscosity_Bender2017::initParameters()
   {
       ViscosityBase::initParameters();
   
       ITERATIONS = createNumericParameter("viscoIterations", "Iterations", &m_iterations);
       setGroup(ITERATIONS, "Fluid Model|Viscosity");
       setDescription(ITERATIONS, "Iterations required by the viscosity solver.");
       getParameter(ITERATIONS)->setReadOnly(true);
   
       MAX_ITERATIONS = createNumericParameter("viscoMaxIter", "Max. iterations (visco)", &m_maxIter);
       setGroup(MAX_ITERATIONS, "Fluid Model|Viscosity");
       setDescription(MAX_ITERATIONS, "Max. iterations of the viscosity solver.");
       static_cast<NumericParameter<unsigned int>*>(getParameter(MAX_ITERATIONS))->setMinValue(1);
   
       MAX_ERROR = createNumericParameter("viscoMaxError", "Max. visco error", &m_maxError);
       setGroup(MAX_ERROR, "Fluid Model|Viscosity");
       setDescription(MAX_ERROR, "Max. error of the viscosity solver.");
       RealParameter* rparam = static_cast<RealParameter*>(getParameter(MAX_ERROR));
       rparam->setMinValue(1e-6);
   }
   
   void Viscosity_Bender2017::step()
   {
       Simulation *sim = Simulation::getCurrent();
       const int numParticles = (int) m_model->numActiveParticles();
       if (numParticles == 0)
           return;
   
       const unsigned int nFluids = sim->numberOfFluidModels();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
       FluidModel *model = m_model;
       const unsigned int fluidModelIndex = model->getPointSetIndex();
       const unsigned int maxIter = m_maxIter;
       const Real maxError = m_maxError;   
       const Real maxError2 = maxError*maxError;
       const Real density0 = model->getDensity0();
   
       const Real h = TimeManager::getCurrent()->getTimeStepSize();
   
       // Compute factors
       computeViscosityFactor();
       computeTargetStrainRate();
   
       m_iterations = 0;
       while (m_iterations < maxIter)
       {
           // Compute viscosity constraint value
           #pragma omp parallel default(shared)
           {
               #pragma omp for schedule(static) nowait 
               for (int i = 0; i < numParticles; i++)
               {
                   const Vector3r &xi = m_model->getPosition(i);
                   const Vector3r &vi = m_model->getVelocity(i);
                   const Real density_i = m_model->getDensity(i);
   
                   Vector6r viscosityC;
                   viscosityC.setZero();
   
                   // Fluid
                   forall_fluid_neighbors_in_same_phase(
                       const Vector3r &vj = m_model->getVelocity(neighborIndex);
                       const Vector3r gradW = sim->gradW(xi - xj);
                       const Vector3r vji = vj - vi;
   
                       const Real m = m_model->getMass(neighborIndex);
                       const Real m2 = m * static_cast<Real>(2.0);
                       viscosityC[0] += m2 * vji[0] * gradW[0];
                       viscosityC[1] += m2 * vji[1] * gradW[1];
                       viscosityC[2] += m2 * vji[2] * gradW[2];
                       viscosityC[3] += m * (vji[0] * gradW[1] + vji[1] * gradW[0]);
                       viscosityC[4] += m * (vji[0] * gradW[2] + vji[2] * gradW[0]);
                       viscosityC[5] += m * (vji[1] * gradW[2] + vji[2] * gradW[1]);
                   )
   
                   viscosityC = (0.5 / density_i) * viscosityC - getTargetStrainRate(i);
    
                   getViscosityLambda(i) = viscosityC;
               }
           }
    
           Real avgStrainRateError = 0.0;
           for (int i = 0; i < (int)numParticles; i++)
               for (unsigned int j = 0; j < 6; j++)
                   avgStrainRateError += fabs(getViscosityLambda(i)[j]);
           avgStrainRateError = (static_cast<Real>(1.0) / (static_cast<Real>(6.0) *(Real)numParticles)) * avgStrainRateError;
   
           // Compute viscosity constraint value
           #pragma omp parallel default(shared)
           {
               #pragma omp for schedule(static) nowait 
               for (int i = 0; i < (int)numParticles; i++)
               {
                   getViscosityLambda(i) = -getViscosityFactor(i) * getViscosityLambda(i);
               }
           }
   
           
           // Apply impulses
           #pragma omp parallel default(shared)
           {
               #pragma omp for schedule(static) nowait 
               for (int i = 0; i < numParticles; i++)
               {
                   const Vector3r &xi = m_model->getPosition(i);
                   Vector3r &vi = m_model->getVelocity(i);
                   const Real density_i = m_model->getDensity(i);
   
   
                   // Fluid
                   Eigen::Matrix<Real, 3, 6> gradT;
                   forall_fluid_neighbors_in_same_phase(
                       const Vector3r gradW = sim->gradW(xi - xj);
                       const Real density_j = m_model->getDensity(neighborIndex);
   
                       gradT.setZero();
                       gradT(0,0) = static_cast<Real>(2.0) * gradW[0];
                       gradT(0,3) = gradW[1];
                       gradT(0,4) = gradW[2];
    
                       gradT(1,1) = static_cast<Real>(2.0) * gradW[1];
                       gradT(1,3) = gradW[0];
                       gradT(1,5) = gradW[2];
    
                       gradT(2,2) = static_cast<Real>(2.0) * gradW[2];
                       gradT(2,4) = gradW[0];
                       gradT(2,5) = gradW[1];
    
                       vi -= m_model->getMass(neighborIndex)* 0.5 * gradT * ((m_model->getMass(neighborIndex) / (density_i*density_i)) * getViscosityLambda(i) +
                           (m_model->getMass(neighborIndex) / (density_j*density_j)) * getViscosityLambda(neighborIndex));
                   )
               }
           }
           m_iterations++;
   
           if (avgStrainRateError < maxError)
               break;
           
       }
       INCREASE_COUNTER("Visco iterations", static_cast<Real>(m_iterations));
   
       // Compute viscosity forces (XSPH) with boundary to simulate simple friction
       const Real invH = (static_cast<Real>(1.0) / h);
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static)  
           for (int i = 0; i < (int)numParticles; i++)
           {
               const Vector3r &xi = m_model->getPosition(i);
               const Vector3r &vi = m_model->getVelocity(i);
               Vector3r &ai = m_model->getAcceleration(i);
               const Real density_i = m_model->getDensity(i);
   
               // Boundary
               if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
               {
                   forall_boundary_neighbors(
                       const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
                       const Vector3r a = -invH * 0.1 * m_viscosity * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi - vj)* sim->W(xi - xj);
                       ai += a;
                       bm_neighbor->addForce(xj, -m_model->getMass(i) * a);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
               {
                   forall_density_maps(
                       Vector3r vj;
                       bm_neighbor->getPointVelocity(xi, vj);
                       const Vector3r a = -invH * 0.1 * m_viscosity * (density0 / density_i) * (vi-vj)* rho;
                       ai += a;
                       bm_neighbor->addForce(xj, -m_model->getMass(i) * a);
                   );
               }
               else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
               {
                   forall_volume_maps(
                       Vector3r vj;
                       bm_neighbor->getPointVelocity(xj, vj);
                       const Vector3r a = -invH * 0.1 * m_viscosity * (density0 * Vj / density_i) * (vi-vj)* sim->W(xi - xj);
                       ai += a;
                       bm_neighbor->addForce(xj, -m_model->getMass(i) * a);
                   );
               }
           }
       }
   }
   
   
   void Viscosity_Bender2017::reset()
   {
   }
   
   
   void Viscosity_Bender2017::computeViscosityFactor()
   {
       // Init parameters
   
       Simulation *sim = Simulation::getCurrent();
       const int numParticles = (int) m_model->numActiveParticles();
       const unsigned int nFluids = sim->numberOfFluidModels();
       const FluidModel *model = m_model;
       const unsigned int fluidModelIndex = model->getPointSetIndex();
   
       #pragma omp parallel default(shared)
       {
           // Compute inverted viscosity matrix
   
           #pragma omp for schedule(static)  
           for (int i = 0; i < numParticles; i++)
           {
               // Compute viscosity matrix
               const Vector3r &xi = m_model->getPosition(i);
               const Real density_i = m_model->getDensity(i);
               Matrix6r &Kinv = getViscosityFactor(i);
               Matrix6r K;
               K.setZero();
       
               Eigen::Matrix<Real, 6, 3> grad_i;
               grad_i.setZero();
   
               // Fluid
               Eigen::Matrix<Real, 6, 3> grad_j;
               forall_fluid_neighbors_in_same_phase(
                   const Vector3r gradW = sim->gradW(xi - xj);
   
                   grad_j.setZero();
                   grad_j(0,0) = static_cast<Real>(2.0) * gradW[0];
                   grad_j(3,0) = gradW[1];
                   grad_j(4,0) = gradW[2];
   
                   grad_j(1,1) = static_cast<Real>(2.0) * gradW[1];
                   grad_j(3,1) = gradW[0];
                   grad_j(5,1) = gradW[2];
   
                   grad_j(2,2) = static_cast<Real>(2.0) * gradW[2];
                   grad_j(4,2) = gradW[0];
                   grad_j(5,2) = gradW[1];
   
                   grad_j = ((0.5 / density_i) * m_model->getMass(neighborIndex)) * grad_j;
                   grad_i -= grad_j;
   
                   Matrix6r Klocal;
                   viscoGradientMultTransposeRightOpt((1.0 / m_model->getDensity(i)) * grad_j, grad_j, Klocal);
                   K += Klocal;
               )
   
   
               Matrix6r Klocal;
               viscoGradientMultTransposeRightOpt((1.0 / m_model->getDensity(i)) * grad_i, grad_i, Klocal);
               K += Klocal;
               
               Vector6r Kdiag_inv;
               for (unsigned l = 0; l < 6; l++)
               {
   
                   if (fabs(K(l,l)) < 1.0e-6)
                       Kdiag_inv[l] = 1.0;
                   else
                       Kdiag_inv[l] = static_cast<Real>(1.0) / K(l,l);
               }
               Matrix6r precondK;
               for (unsigned k = 0; k < 6; k++)
               {
                   for (unsigned l = 0; l < 6; l++)
                   {
                       precondK(k,l) = Kdiag_inv[k] * K(k,l);
                   }
               }
   
               if (fabs(precondK.determinant()) < 1.0e-6)
                   Kinv.setZero();
               else
                   Kinv = precondK.inverse();
   
               for (unsigned k = 0; k < 6; k++)
               {
                   for (unsigned l = 0; l < 6; l++)
                   {
                       Kinv(k,l) = Kinv(k,l) * Kdiag_inv[l];
                   }
               }
           }
       }
   }
   
   void Viscosity_Bender2017::computeTargetStrainRate()
   {
       Simulation *sim = Simulation::getCurrent();
       const int numParticles = (int) m_model->numActiveParticles();
       const unsigned int fluidModelIndex = m_model->getPointSetIndex();
       const unsigned int nFluids = sim->numberOfFluidModels();
       FluidModel *model = m_model;
           
       // Compute target strain rate
       #pragma omp parallel default(shared)
       {
           #pragma omp for schedule(static) nowait 
           for (int i = 0; i < numParticles; i++)
           {
               const Vector3r &xi = m_model->getPosition(i);
               const Vector3r &vi = m_model->getVelocity(i);
               const Real density_i = m_model->getDensity(i);
   
               Vector6r &strainRate = getTargetStrainRate(i);
               strainRate.setZero();
   
               // Fluid
               forall_fluid_neighbors_in_same_phase(
                   const Vector3r &vj = m_model->getVelocity(neighborIndex);
   
                   const Vector3r gradW = sim->gradW(xi - xj);
                   const Vector3r vji = vj - vi;
                   const Real m = m_model->getMass(neighborIndex);
                   const Real m2 = m * static_cast<Real>(2.0);
                   strainRate[0] += (static_cast<Real>(1.0)-m_viscosity) * m2 * vji[0] * gradW[0];
                   strainRate[1] += (static_cast<Real>(1.0)-m_viscosity) * m2 * vji[1] * gradW[1];
                   strainRate[2] += (static_cast<Real>(1.0)-m_viscosity) * m2 * vji[2] * gradW[2];
                   strainRate[3] += (static_cast<Real>(1.0)-m_viscosity) * m * (vji[0] * gradW[1] + vji[1] * gradW[0]);
                   strainRate[4] += (static_cast<Real>(1.0)-m_viscosity) * m * (vji[0] * gradW[2] + vji[2] * gradW[0]);
                   strainRate[5] += (static_cast<Real>(1.0)-m_viscosity) * m * (vji[1] * gradW[2] + vji[2] * gradW[1]);
               )
               strainRate = (static_cast<Real>(0.5) / density_i) * strainRate;
           }
       }
   }
   
   void Viscosity_Bender2017::performNeighborhoodSearchSort()
   {
       const unsigned int numPart = m_model->numActiveParticles();
       if (numPart == 0)
           return;
   
       Simulation *sim = Simulation::getCurrent();
       auto const& d = sim->getNeighborhoodSearch()->point_set(m_model->getPointSetIndex());
       d.sort_field(&m_targetStrainRate[0]);
       d.sort_field(&m_viscosityFactor[0]);
       d.sort_field(&m_viscosityLambda[0]);
   }