testKernelGrav.c 3.36 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *                    James Willis (james.s.willis@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
20
#include "const.h"
21
#include "kernel_gravity.h"
22
#include "kernel_long_gravity.h"
23 24 25 26 27 28

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>

29
#define numPoints (1 << 7)
30 31 32 33

/**
 * @brief The Gadget-2 gravity kernel function
 *
34 35
 * Taken from Gadget-2.0.7's forcetree.c lines 2755-2800
 *
36
 * @param r The distance between particles
37
 * @param epsilon The cut-off distance of the kernel
38
 */
39 40 41 42 43 44 45 46 47 48 49 50
float gadget(float r, float epsilon) {

  const float h = epsilon;
  const float h_inv = 1.f / h;

  const float u = r * h_inv;

  if (u >= 1) {
    const float r_inv = 1. / r;

    return r_inv * r_inv * r_inv;
  } else {
51
    if (u < 0.5)
52 53
      return h_inv * h_inv * h_inv *
             (10.666666666667 + u * u * (32.0 * u - 38.4));
54
    else
55 56 57
      return h_inv * h_inv * h_inv *
             (21.333333333333 - 48.0 * u + 38.4 * u * u -
              10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
58 59 60
  }
}

61
int main(int argc, char *argv[]) {
62 63

  const float h = 3.f;
64
  const float r_max = 6.f;
65 66 67 68 69 70

  for (int k = 1; k < numPoints; ++k) {

    const float r = (r_max * k) / numPoints;
    const float gadget_w = gadget(r, h);

71 72 73 74
    const float h_inv = 1.f / h;
    const float h_inv3 = h_inv * h_inv * h_inv;
    const float u = r * h_inv;

75
    float swift_w;
76
    if (r >= h) {
77
      swift_w = 1 / (r * r * r);
78 79 80 81

    } else {
      kernel_grav_eval(u, &swift_w);
      swift_w *= h_inv3;
82 83
    }

84
    if (fabsf(gadget_w - swift_w) > 1e-5 * fabsf(gadget_w)) {
85 86

      printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
Matthieu Schaller's avatar
Matthieu Schaller committed
87
             gadget_w, swift_w);
88

89 90 91 92 93 94
      printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
      return 1;
    }
  }

  printf("\nAll values are consistent\n");
95 96

  /* Now test the long range function */
97
  /* const float a_smooth = 4.5f; */
98

99
  /* for (int k = 1; k < numPoints; ++k) { */
100

101
  /*   const float r = (r_max * k) / numPoints; */
102

103
  /*   const float u = r / a_smooth; */
104

105 106
  /*   float swift_w; */
  /*   kernel_long_grav_eval(u, &swift_w); */
107

108
  /*   float gadget_w = erfcf(u / 2) + u * expf(-u * u / 4) / sqrtf(M_PI); */
109

110
  /*   if (fabsf(gadget_w - swift_w) > 1e-4 * fabsf(gadget_w)) { */
111

112 113 114
  /*     printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r,
   * a_smooth, */
  /*            u, swift_w, gadget_w); */
115

116 117 118 119 120
  /*     printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
   */
  /*     return 1; */
  /*   } */
  /* } */
121

122 123
  return 0;
}