volume.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_ENERGY_VOLUME_HPP
#define PECLET_VORO_ENERGY_VOLUME_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/voro/energy/route.hpp"
#include "peclet/voro/tessellation_view.hpp"
 
namespace peclet::voro::energy {
 
template <class Real>
KOKKOS_INLINE_FUNCTION Real dTargetVolume(Real V, Real Vref) {
  return Real(2) * (V / Vref - Real(1)) / Vref;
}
template <class Real>
KOKKOS_INLINE_FUNCTION Real dLogBarrier(Real V, Real mu) {
  return V > Real(0) ? -mu / V : Real(0);
}
template <class Real>
KOKKOS_INLINE_FUNCTION Real dFreeEnergy(Real V, Real Vref) {
  return V > Real(0) ? -Vref / V : Real(0);
}
 
template <class Real, class Policy = Voronoi, class Sdf = NoSdf>
void volumeGradientForce(const TessellationView<Real>& view,
                         const Kokkos::View<Real*, peclet::core::MemSpace>& dEdV,
                         const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                         const Kokkos::View<Real*, peclet::core::MemSpace>& weight, Real L,
                         const Kokkos::View<Real*, peclet::core::MemSpace>& force,
                         const Kokkos::View<Real*, peclet::core::MemSpace>& forceW = {},
                         const Sdf& sdf = Sdf{}) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Kokkos::parallel_for(
      "energy.volume", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        const Real de = dEdV(i);
        if (de == Real(0))
          return;
        Real gl[detail::kMaxLocalFacets * 3] = {};  // value-init: -Wmaybe-uninitialized
        const int b = view.facetBegin(i), nf = view.facetEnd(i) - b;
        if (nf > detail::kMaxLocalFacets)
          return;
        const bool fd = view.hasWallFD();  // exact wall part (sdfWallFD) instead of the chain
        for (int k = 0; k < nf; ++k) {
          const int f = b + k;
          if (view.facetNbr(f) == kBoundaryFacet) {
            for (int c = 0; c < 3; ++c)
              gl[3 * k + c] = fd ? Real(0) : de * view.area(f, c);
          } else {
            for (int c = 0; c < 3; ++c)
              gl[3 * k + c] = Real(2) * de * view.connect(f, c);
          }
        }
        detail::routeCellGradient<Real, Policy>(view, i, gl, nf, pos, weight, L, sdf, force,
                                                forceW);
        if (fd)
          for (int c = 0; c < 3; ++c)
            Kokkos::atomic_add(&force(3 * i + c), de * view.wallDV(i, c));
      });
}
 
}  // namespace peclet::voro::energy
 
#endif  // PECLET_VORO_ENERGY_VOLUME_HPP