SPHKernels.h

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <array>
#include <cmath>
 
#include "autopas/utils/ArrayMath.h"
 
namespace sphLib {
class SPHKernels {
  static constexpr double pi = M_PI;  // atan(1.0) * 4.0;
  static constexpr double kernelSupportRadius = 2.5;
  // const double C_CFL = 0.3;
 public:
  static inline double getKernelSupportRadius() { return kernelSupportRadius; }
 
  static inline double W(const std::array<double, 3> &dr, const double h) {
    const double dr2 = autopas::utils::ArrayMath::dot(dr, dr);
    return W(dr2, h);
  }
 
  static inline double W(const double dr2, const double h) {
    const double H = kernelSupportRadius * h;
    if (dr2 < H * H) {
      const double s = sqrt(dr2) / H;  // sqrt(dr * dr) / H;
      const double s1 = 1.0 - s;       // (1.0 - s < 0) ? 0 : 1.0 - s;
      const double s2 = std::max(0., 0.5 - s);
      double r_value = (s1 * s1 * s1) - 4.0 * (s2 * s2 * s2);  //  pow(s1, 3) - 4.0 * pow(s2, 3);
      // if # of dimension == 3
      r_value *= 16.0 / pi / (H * H * H);
      return r_value;
    } else {
      return 0.;
    }
  }
 
  static unsigned long getFlopsW();
 
  static inline std::array<double, 3> gradW(const std::array<double, 3> &dr, const double h) {
    using namespace autopas::utils::ArrayMath::literals;
 
    const double H = kernelSupportRadius * h;
    const double drabs = autopas::utils::ArrayMath::L2Norm(dr);
    const double s = drabs / H;  // sqrt(dr * dr) / H;
    const double s1 = (1.0 - s < 0) ? 0 : 1.0 - s;
    const double s2 = (0.5 - s < 0) ? 0 : 0.5 - s;
    double r_value = -3.0 * (s1 * s1) + 12.0 * (s2 * s2);  // - 3.0 * pow(s1, 2) + 12.0 * pow(s2, 2);
    // if # of dimension == 3
    r_value *= 16.0 / pi / (H * H * H);
    const double scale = r_value / (drabs * H + 1.0e-6 * h);
    return dr * scale;  // dr * r_value / (sqrt(dr * dr) * H + 1.0e-6 * h);
  }
};
}  // namespace sphLib