viscous.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_PHYSICS_VISCOUS_HPP
#define PECLET_VORO_PHYSICS_VISCOUS_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/tessellation_view.hpp"
 
namespace peclet::voro {
namespace physics {
 
template <class Real>
void velocityGradient(const TessellationView<Real>& view,
                      const Kokkos::View<int*, peclet::core::MemSpace>& recip,
                      const Kokkos::View<int*, peclet::core::MemSpace>& cellOfFacet,
                      const Kokkos::View<Real*, peclet::core::MemSpace>& vel,
                      const Kokkos::View<Real*, peclet::core::MemSpace>& grad) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Kokkos::parallel_for(
      "visc.grad", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        Real g[9];
        for (int t = 0; t < 9; ++t)
          g[t] = 0;
        const Real vi[3] = {vel(3 * i + 0), vel(3 * i + 1), vel(3 * i + 2)};
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f) {
          const int r = recip(f);
          if (r < 0)
            continue;  // boundary facet
          const int nb = cellOfFacet(r);
          for (int p = 0; p < 3; ++p) {
            const Real vn = vel(3 * nb + p);
            for (int k = 0; k < 3; ++k)
              g[3 * p + k] += vi[p] * view.connect(f, k) - vn * view.connect(r, k);
          }
        }
        const Real invV = Real(1) / view.volume(i);
        for (int t = 0; t < 9; ++t)
          grad(9 * i + t) = g[t] * invV;
      });
}
 
template <class Real>
void viscousStress(const Kokkos::View<Real*, peclet::core::MemSpace>& grad,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& visc,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& bulkVisc,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& stress) {
  using Exec = peclet::core::ExecSpace;
  const int N = static_cast<int>(visc.extent(0));
  Kokkos::parallel_for(
      "visc.stress", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        const Real mu = visc(i), mub = bulkVisc(i);
        Real gk[9];
        for (int t = 0; t < 9; ++t)
          gk[t] = grad(9 * i + t);
        Real divv = gk[0] + gk[4] + gk[8];
        Real s[9];
        for (int k = 0; k < 3; ++k)
          for (int m = 0; m < 3; ++m)
            s[3 * k + m] = mu * (gk[3 * k + m] + gk[3 * m + k]);
        const Real lam = (mub - Real(2) / Real(3) * mu) * divv;
        s[0] += lam;
        s[4] += lam;
        s[8] += lam;
        for (int t = 0; t < 9; ++t)
          stress(9 * i + t) = s[t];
      });
}
 
template <class Real>
void addViscousForce(const TessellationView<Real>& view,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& stress,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& force) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Kokkos::parallel_for(
      "visc.force", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        Real f[3] = {0, 0, 0};
        for (int ff = view.facetBegin(i); ff < view.facetEnd(i); ++ff) {
          const int nb = view.facetNbr(ff);
          if (nb < 0)
            continue;  // boundary facet
          for (int k = 0; k < 3; ++k)
            for (int m = 0; m < 3; ++m)
              f[k] +=
                  (stress(9 * nb + 3 * k + m) - stress(9 * i + 3 * k + m)) * view.connect(ff, m);
        }
        for (int k = 0; k < 3; ++k)
          force(3 * i + k) += f[k];
      });
}
 
template <class Real>
void viscousForce(const TessellationView<Real>& view,
                  const Kokkos::View<int*, peclet::core::MemSpace>& recip,
                  const Kokkos::View<int*, peclet::core::MemSpace>& cellOfFacet,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& vel,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& visc,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& bulkVisc,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& force,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& grad,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& stress) {
  velocityGradient(view, recip, cellOfFacet, vel, grad);
  viscousStress(grad, visc, bulkVisc, stress);
  addViscousForce(view, stress, force);
}
 
template <class Real>
void viscousForce(const TessellationView<Real>& view,
                  const Kokkos::View<int*, peclet::core::MemSpace>& recip,
                  const Kokkos::View<int*, peclet::core::MemSpace>& cellOfFacet,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& vel,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& visc,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& bulkVisc,
                  const Kokkos::View<Real*, peclet::core::MemSpace>& force) {
  const int N = view.numCells();
  Kokkos::View<Real*, peclet::core::MemSpace> grad("visc.grad", (size_t)9 * N);
  Kokkos::View<Real*, peclet::core::MemSpace> stress("visc.stress", (size_t)9 * N);
  viscousForce(view, recip, cellOfFacet, vel, visc, bulkVisc, force, grad, stress);
}
 
}  // namespace physics
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_PHYSICS_VISCOUS_HPP