kernel_hydro.h

Go to the documentation of this file.

#pragma once
#include "cmath"
#include "tarch/tarch.h"
 
#include "peano4/utils/Globals.h"
#include "tarch/compiler/CompilerSpecificSettings.h"
 
 
namespace exaparty {
  namespace kernels {
    namespace legacy {
      namespace kernelHydro {
 
#if defined(QUARTIC_SPLINE_KERNEL)
        /* -------------------------------------------------------------------------
         */
 
        /* Quartic spline M5 kernel -- SWIFT version */
 
#define kernel_name "Quartic spline (M5)"
/* Coefficients for the kernel. */
#define kernel_degree 4 /* Degree of the polynomial */
#define kernel_ivals 5  /* Number of branches */
 
#if HYDRO_DIMENSION == 1
#define kernel_gamma ((float)(1.936492))
#define kernel_constant ((float)(3125. / 768.))
#elif HYDRO_DIMENSION == 2
#define kernel_gamma ((float)(1.977173))
#define kernel_constant ((float)(46875. * M_1_PI / 2398.))
#elif HYDRO_DIMENSION == 3
#define kernel_gamma ((float)(2.018932))
#define kernel_constant ((float)(15625. * M_1_PI / 512.))
#else
#pragma error "Hydro dimension undefined"
#endif
 
#if defined(GPUOffloadingOMP)
#pragma omp begin declare target
#endif
        static const double
          kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
          __attribute__((aligned(AlignmentOnHeap)))
          = {
            6.0,  0.0,  -2.4, 0.0,  0.368, /* 0 < u < 0.2 */
            -4.0, 8.0,  -4.8, 0.32, 0.352, /* 0.2 < u < 0.4 */
            -4.0, 8.0,  -4.8, 0.32, 0.352, /* 0.4 < u < 0.6 */
            1.0,  -4.0, 6.0,  -4.0, 1.0,   /* 0.6 < u < 0.8 */
            1.0,  -4.0, 6.0,  -4.0, 1.0,   /* 0.8 < u < 1 */
            0.0,  0.0,  0.0,  0.0,  0.0};  /* 1 < u */
#if defined(GPUOffloadingOMP)
#pragma omp end declare target
#endif
 
/* ------------------------------------------------------------------------- */
 
// @todo implement additional kernels (M4)
#else
 
#error "An SPH projection kernel function must be chosen !"
/* ------------------------------------------------------------------------- */
#endif
 
/* ------------------------------------------------------------------------- */
 
/* Define generic values (i.e. common to all kernels) */
 
/* First some powers of gamma = H/h */
#define kernel_gamma_inv ((float)(1. / kernel_gamma))
#define kernel_gamma2 ((float)(kernel_gamma * kernel_gamma))
 
/* define gamma^d, gamma^(d+1), 1/gamma^d and 1/gamma^(d+1) */
#if HYDRO_DIMENSION == 1
#define kernel_gamma_dim ((float)(kernel_gamma))
#define kernel_gamma_dim_plus_one ((float)(kernel_gamma * kernel_gamma))
#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma)))
#define kernel_gamma_inv_dim_plus_one \
  ((float)(1. / (kernel_gamma * kernel_gamma)))
#elif HYDRO_DIMENSION == 2
#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma))
#define kernel_gamma_dim_plus_one \
  ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma * kernel_gamma)))
#define kernel_gamma_inv_dim_plus_one \
  ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
#elif HYDRO_DIMENSION == 3
#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
#define kernel_gamma_dim_plus_one \
  ((float)(kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma))
#define kernel_gamma_inv_dim \
  ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
#define kernel_gamma_inv_dim_plus_one \
  ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma)))
#endif
 
/* The number of branches (floating point conversion) */
#define kernel_ivals_f ((float)(kernel_ivals))
 
/* Kernel self contribution (i.e. W(0,h)) */
#define kernel_root \
  ((float)(kernel_coeffs[kernel_degree]) * kernel_constant \
   * kernel_gamma_inv_dim)
 
        /* -------------------------------------------------------------------------
         */
 
        static void kernel_eval(double u, double& W) {
 
          /* Go to the range [0,1[ from [0,H[ */
          const double x = u * kernel_gamma_inv;
 
          /* Pick the correct branch of the kernel */
          const int temp = (int)(x * kernel_ivals_f);
          const int ind  = temp > kernel_ivals ? kernel_ivals : temp;
          const double* const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
 
          /* First two terms of the polynomial ... */
          double w = coeffs[0] * x + coeffs[1];
 
          /* ... and the rest of them */
          for (int k = 2; k <= kernel_degree; k++)
            w = x * w + coeffs[k];
 
          w = std::max(w, 0.0);
 
          /* Return everything */
          W = w * kernel_constant * kernel_gamma_inv_dim;
        }
 
        static void kernel_deval(double u, double& W, double& dW_dx) {
          /* Go to the range [0,1[ from [0,H[ */
          const double x = u * kernel_gamma_inv;
 
          /* Pick the correct branch of the kernel */
          const int temp = (int)(x * kernel_ivals_f);
          const int ind  = temp > kernel_ivals ? kernel_ivals : temp;
          const double* const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
 
          /* w term in Horner's scheme */
          double w = coeffs[0];
          for (int k = 1; k <= kernel_degree; k++) {
            w = w * x + coeffs[k];
          }
 
          /* dw_dx term in Horner's scheme */
 
          double dw_dx = coeffs[0] * kernel_degree;
 
          /* ... and the rest of them */
          for (int k = 1; k <= kernel_degree - 1; k++) {
            dw_dx = dw_dx * x + coeffs[k] * (kernel_degree - k);
          }
 
          w     = std::max(w, 0.0);
          dw_dx = std::min(dw_dx, 0.0);
 
          /* Return everything */
          W     = u == 0 ? 0: w * kernel_constant * kernel_gamma_inv_dim;
          dW_dx = u == 0 ? 0: dw_dx * kernel_constant * kernel_gamma_inv_dim_plus_one;
        }
      } // namespace kernelHydro
    }   // namespace legacy
  }     // namespace kernels
} // namespace exaparty