euler_pressure.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_PHYSICS_EULER_PRESSURE_HPP
#define PECLET_VORO_PHYSICS_EULER_PRESSURE_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 eulerPressureForce(const TessellationView<Real>& view,
                        const Kokkos::View<int*, peclet::core::MemSpace>& recip,
                        const Kokkos::View<int*, peclet::core::MemSpace>& cellOfFacet, Real pressEq,
                        Real volAvg, const Kokkos::View<Real*, peclet::core::MemSpace>& force) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Kokkos::parallel_for(
      "euler.pressure", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        const Real pressI = pressEq * volAvg / view.volume(i);
        Real fx = 0, fy = 0, fz = 0;
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f) {
          // self half-contribution (always present, including boundary facets)
          fx -= pressI * view.connect(f, 0);
          fy -= pressI * view.connect(f, 1);
          fz -= pressI * view.connect(f, 2);
          const int g = recip(f);
          if (g < 0)
            continue;  // boundary / no reciprocal
          const Real pressJ = pressEq * volAvg / view.volume(cellOfFacet(g));
          fx += pressJ * view.connect(g, 0);
          fy += pressJ * view.connect(g, 1);
          fz += pressJ * view.connect(g, 2);
        }
        force(3 * i + 0) = fx;
        force(3 * i + 1) = fy;
        force(3 * i + 2) = fz;
      });
}
 
}  // namespace physics
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_PHYSICS_EULER_PRESSURE_HPP